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OBJECTIVE 


The  goals  of  this  effort  are: 

•  Select  and  develop  methods  for  the  analysis  of  the  thermomechanical  behavior  of  tex¬ 
tile  composite  materials. 

•  Perform  analysis  of  selected  textile  forms,  including  2-D  braids  and  weaves,  and  poten¬ 
tially  3-D  architectures. 

•  Develop  failure  initiation  and  growth  methodology  for  textile  composites,  including  the 
growth  and  coalescence  of  internal  cracking. 

•  Study  the  relationship  of  fabric  architecture  and  failure  sequence  in  a  variety  of  textile 
composites. 

SUMMARY  OF  ACCOMPLISHMENTS 

To  date,  different  analytical  techniques  were  examined  for  their  suitability  for  thermome¬ 
chanical  analysis  of  textile  composites. 

•  The  concept  of  equivalent  homogeneity  was  used  to  derive  relationships  between  av¬ 
erage  stresses  and  average  strains  within  a  fabric  (see  Appendix  C).  Representative 
volume  elements  (RVE)  consisting  of  several  non-homogeneous  layers  (fabric/matrix) 
were  considered.  Methods  for  calculating  simple  bounds  on  the  effective  properties  of 
woven  laminates  were  also  suggested. 

•  Various  numerical  approaches  for  woven  composites  using  h-element,  p-element  and 
novel  meshing  techniques  were  investigated  (see  Appendix  D).  Automated  mesh  gen¬ 
eration  codes  were  developed  to  rapidly  create  detailed  FE  models  of  typical  two- 
dimensional  textile  forms.  The  mesh  generation  tools  are  modular  and  allow  parametric 
modeling  of  the  fabric  geometry.  The  symmetry  and  antisymmetry  of  material  distribu¬ 
tion  arising  out  of  the  repeating  weave  pattern  can  be  used  advantageously  to  reduce 
the  problem  size.  This  has  been  demonstrated  in  the  case  of  plain  woven  fabrics.  The 
models  were  analyzed  using  a  commercial  FEA  code,  ABAQUS,  to  compute  their  elas¬ 
tic  properties  and  their  mechanical  response. 
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•  Nonlinear  material  behavior  and  damage  growth  in  textiie  composites  was  investigated 
using  damage  mechanics  based  constitutive  laws.  The  strain  energy  dissipation  was 
chosen  as  the  damage  parameter  and  used  to  describe  a  set  of  damage  surfaces. 
Materiai  models  were  developed  for  each  of  the  weave  constituents  viz.  yarn  and  pure 
matrix  to  describe  their  constitutive  relationships  and  damage  state.  Interface  subrou¬ 
tines  were  deveioped  to  iink  the  deveioped  material  model  codes  with  a  generai- 
purpose  finite  element  analysis  program. 

•  The  developed  damage  modeling  methodology  is  demonstrated  for  plain  woven  fabrics. 
Progressive  failure  analyses  of  some  plain  weave  composites  subjected  to  tension  and 
in-plane  shear  are  carried  out.  The  different  damage  mechanisms  causing  failure  of  the 
fabric  under  on-axes  and  off-axes  loading  conditions  are  studied.  Parametric  studies  to 
study  the  effect  of  yam  waviness  and  fiber  volume  fraction  on  the  predicted  elastic 
properties  of  plain  weaves  are  also  done. 
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INTRODUCTION 


Micromechanics  of  textile  composites  has  been  the  focus  of  study  by  many  investigators  [1- 
7].  A  survey  of  the  literature  in  this  area  is  documented  in  Ref.  [8].  The  intricate  geometry  of  the 
yarns  and  varying  spatial  material  properties  within  the  fabric  makes  the  response  of  such 
composites  different  from  that  of  unidirectional  composites.  Most  of  the  analytical  studies  on 
such  composites  have  been  focussed  on  predicting  its  elastic  properties  [1-7].  Studies  relating 
to  prediction  of  their  strength  and  failure  modes  are  relatively  few  [9-12].  To  understand  the  be¬ 
havior  and  damage  mechanisms  causing  failure  of  fabric  composites,  a  detailed  modeling  of 
the  geometry  and  material  properties  is  necessary.  Such  modeling  is  possible  by  the  use  of 
three-dimensional  finite  elements.  Existing  modeling  tools  and  FE  preprocessors  are  not  ade¬ 
quate  for  rapidly  creating  complex  3-dimensional  models  needed  for  the  purpose.  Therefore, 
mesh  generation  programs  were  developed  which  would  easily  interface  with  a  general  pur¬ 
pose  FE  analysis  software,  ABAQUS  [13].  FE  models  can  be  generated  for  plain  weaves  and 
higher  order  satin-weaves  using  a  minimal  definition  of  the  textile  architecture. 

Nonlinear  material  behavior  of  composites  is  due  to  damage  accumulation,  which  causes 
changes  in  the  stiffnesses  of  the  material.  Material  nonlinearity  combined  with  high  local  stress 
fields  can  result  in  a  variety  of  damage  mechanisms  such  as  matrix  cracking,  fiber  breakage, 
and  delamination  which  influence  the  overall  properties  of  the  composite.  It  is  well  known  that 
macroscopic  failure  is  preceded  by  an  accumulation  of  different  types  of  microscopic  damage 
and  stiffness  losses  due  to  the  accumulation  of  such  damage,  which  cause  significant  load  re¬ 
distribution.  Final  failure  usually  occurs  due  to  the  development  of  one  or  more  areas  of  highly 
localized  stiffness  losses,  which  behave  as  macroscopic  damages  or  cracks.  Hence,  a  proper 
modeling  of  the  damage  growth,  stiffness  loss  and  load  redistribution  is  essential  for  a  reliable 
prediction  of  the  weave  response. 

In  the  present  approach,  damages  in  the  composite  constituents  are  modeled  on  a  contin¬ 
uum  basis  and  related  to  its  material  constitutive  behavior.  The  3D  constitutive  laws  describing 
matrix  and  yarn  behavior  are  developed  using  a  damage  mechanics  based  approach  with  the 
dissipated  energy  density  as  the  damage  parameter.  The  strain  energy  dissipation  (SED)  con¬ 
cept  is  employed  to  describe  the  damage  state  and  current  stiffnesses  of  the  weave  constitu¬ 
ents.  The  yarns  are  treated  to  be  transversely  isotropic  while  the  pure  matrix  pockets  are  as¬ 
sumed  to  be  isotropic. 
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The  developed  methodology  for  geometric  and  damage  modeling  in  fabric  composites  is 
demonstrated  through  an  investigation  of  progressive  failure  of  plain  woven  composites  sub¬ 
jected  to  external  loads.  Two  loading  conditions  are  considered  viz.  in-plane  tension  and  in- 
plane  shear.  For  a  realistic  modeling  of  the  weave  response,  both  geometric  and  material  non- 
linearities  are  considered.  The  inititiation  and  growth  of  damage  within  the  fabric  is  investigated 
and  the  important  stresses  causing  failure  identified. 

FABRIC  GEOMETRY 

For  a  realistic  prediction  of  the  fabric  response,  a  proper  modeling  of  the  fabric  architecture 
is  essential.  Within  a  woven  fabric,  the  interlaced  nature  of  the  yarns  causes  an  asymmetry  in 
distribution  of  the  warp  and  fill.  On  one  side  of  the  woven  fabric,  there  is  a  preponderance  of 
warp  while  the  same  is  true  for  the  fill  on  the  other  side.  This  asymmetry  in  material  distribution 
causes  a  coupling  of  the  in-plane  and  out-of-plane  actions.  The  undulations  of  the  yarns  create 
pockets  of  space  within  the  fabric,  which  are  filled  by  the  binding  matrix.  In  the  current  analysis, 
the  yarns  are  assumed  to  follow  a  sinusoidal  path,  over  yarns  running  orthogonally  to  them. 
The  yarn  path  and  cross-sectional  geometry  of  the  yarns  are  similar  to  that  described  in  Ref. 
[3].  The  cross-section  profile  of  the  yarns  and  the  undulation  path  are  described  in  terms  of 
yarn  width,  yarn  count  and  fiber  volume  fraction.  They  provide  an  adequate  representation  of 
the  varying  cross-sectional  profile  and  undulation  of  the  yarns. 

The  periodic  nature  of  the  yarn  undulations  makes  it  possible  to  identify  repetitive  portion  of 
the  woven  fabric.  The  mechanical  response  of  the  fabric  can  then  be  described  by  analyzing 
such  repetitive  portions,  termed  as  unit  cell,  of  the  textile.  The  unit  cell  with  appropriate  bound¬ 
ary  conditions  can  thus  be  used  to  examine  the  behavior  of  the  entire  textile  subjected  to  vari¬ 
ous  loading  conditions.  Unit  cells  of  the  textile  forms  considered  in  the  present  study  are  shown 
in  Figure  1 . 
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(c) 

Figure  1 .  (a)  Plain  Weave  (one  quarter  highlighted)  (b)  5  Harness  Satin  (c)  8  Harness  Satin 
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Finite  Element  Discretization 


The  weave  is  analyzed  using  ABAQUS  [13]  employing  20  noded,  15  noded  and  10  noded 
three-dimensional  elements.  An  in-house  parametric  finite  element  mesh  generator  described 
earlier  permits  an  automatic  meshing  of  the  weave  unit  cell  and  ABAQUS  input  file  generation. 

To  simulate  different  loading  conditions,  appropriate  boundary  conditions  are  used  (see 
Appendix  D).  For  the  case  of  plain  weaves,  symmetry  within  the  periodic  cell  (Figure  2)  can  be 
used  to  identify  the  repeating  unit  cell  (Figure  3).  Further,  the  point  symmetrical  material  distri¬ 
bution  about  the  z-axis  within  the  unit  cell  (Figure  3)  makes  it  possible  to  analyze  only  one 
quarter  (Figure  4)  with  appropriate  boundary  conditions.  This  results  in  a  reduction  of  problem 
size  and  computational  time. 


Figure  2.  Plain  Weave  Periodic  Cell  Figure  3.  Plain  Weave  Unit  Cell 

(some  matrix  portions  (one  quarter  highlighted) 

removed  for  clarity) 


Figure  4.  Quarter  of  Plain  Weave  Unit  Cell 
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Details  of  the  boundary  conditions  applicable  to  the  quarter  cell  plain  weave  model  (Figure 
5)  are  described  in  Appendix  A  and  are  summarized  below; 


Figure  5.  Quarter  Cell  of  Plain  Weave 
For  the  case  of  in-plane  extenslonal  loading  in  the  x  direction: 


Ontheplanex  =  0  : 


Along  they- axis  (0,y,0)  : 
Ontheplaney  =  0  : 


Along  the  X- axis  (x,0,0)  : 
Along  the  z  -  axis  (0,0,  z) : 


On  the  plane  x  =  b/2  : 
On  the  plane  y  =  b  /  2: 


u(0,y,z)  =  -u(0,y,-z) 
v(0,y,z)  =  v(0,y,-z) 
w(0,y,z)  =  -w(0,y,-z) 

u(0,y,0)  =  0.0 

u(x,0,z)  =  u(x,0,-z) 
v(x,0,z)  =  -v(x,0,-z) 
w(x,0,  z)  =  -w(x,0,-z) 

v(x,0,0)  =  0.0 

u(0,0,z)  =  0 
v(0,0,z)  =  0 
w(0,0,z)  =  -w(0,0,-z) 


/h  X  uo 

u(-,y,^)  =  y 


v(x,— ,z)  =  uniform 


(1) 
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The  boundary  conditions  for  in-plane  shear  ioading  are: 


Ontheplanex  =  0 


Ontheplaney  =  0 


On  the  plane  y=— 


Ontheplanex=— 


u(0,y,z)  =  u(0,y,-z) 
v(0,y,z)  =  -v(0,y,-z) 
w(0,y,z)=  w(0,y,-z) 
u(x,0,z)  =  -u(x,0,-z) 
v(x,0,z)  =  v(x,0,-z) 
w(x,0,z)  =  w(x,0,-z) 

w(x,|-,z)  =  0 

v(fy,x)  =  f 


w(-,y,z)  =  0 


(2) 


ELASTIC  ANALYSIS 

The  boundary  conditions  shown  in  Eqns.  1-2,  aiong  with  those  described  eariier  (Appen¬ 
dix  D),  were  used  to  compute  the  elastic  properties  of  plain,  5-harness  satin  and  8-harness 
weaves  to  vaiidate  the  FE  model.  For  the  sake  of  comparison  with  published  analytical  and  ex¬ 
perimental  results  [3],  AS4  Graphite/Epoxy  6501  woven  fabrics  with  a  yarn  fiber  volume  fraction 
of  0.75  were  chosen.  The  width  of  the  yarns  was  1.411  mm  and  its  thickness  was  0.01  mm. 
The  filament  count  in  both  the  warp  and  fiii  tows  was  assumed  to  be  3000  [3].  The  material 
properties  of  AS4-Graphite  fiber  and  epoxy  matrix  shown  in  Tabie  1  were  used  to  compute  the 
effective  elastic  properties  of  the  individual  yarns  using  NDPROP  [14]  and  are  listed  in  Table  2. 
These  effective  properties  of  the  yarns  were  used  to  compute  the  textile  elastic  stiffnesses.  Ta¬ 
ble  3  shows  a  comparison  of  the  current  results  with  other  anaiyticai  and  experimental  results 
for  plain,  5-harness  and  8-harness  satin  weaves.  For  an  AS4  graphite/epoxy  cross-piy  laminate 
having  an  overall  fiber  volume  fraction  equai  to  that  of  the  aforementioned  weaves  (Vf  =  0.64), 
the  elastic  modulus  computed  using  the  CLT  theory  is  67.60  GPa  [3].  It  can  be  seen  that  the 
undulation  of  the  yarns  causes  a  reduction  in  the  in-piane  extensional  stiffness  of  the  weaves 
and  this  reduction  is  maximum  in  plain  weaves  where  the  effect  of  yarn  cross-over  is  maximum. 
The  effect  of  mesh  discretization  on  the  predicted  properties  was  studied  in  the  case  of  plain 
weaves  as  shown  in  Tabie  4. 
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Table  1.  Material  Properties  of  Fiber  and  Matrix  [3] 


Tabie  2.  Materiai  Properties  of  the  Yarns 


Material 

En 

GPa 

E22 

GPa 

V12 

V23 

Gi2 

GPa 

G23 

GPa 

AS4/Epoxy  6501 
(Vf=0.64) 

144.8 

11.73 

0.230 

0.30 

5.52 

3.3 

Table  3.  Comparison  of  Present  Resuits  with  Anaiytical  and  Experimental  Studies 


AS4/Epoxy 

Plain  weave 
(Vf°  =  0.64) 

Present 

Experimental 

[3] 

FEM  [3] 

AS4/Epoxy 

5-harness 

satin 

(Vt°  =  0.64 ) 

Present 

FEM  [3] 

AS4  /  Epoxy 

8-hamess 

satin 

(Vf  =  0.64 ) 

Present 

FEM  [3] 

E. 

GPa  (Msi) 


11.84 

(172) 


0.060  0.420 
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Table  4.  Mesh  Refinement  Studies  (AS4/Epoxy  3501-6  Plain  Weave) 


Number  of 
elements 

Ex 

MPa 

Periodic  cell 

400 

62.338 

28.930 

1296 

62.117 

273.07 

2704 

61.910 

2985.5 

Quarter  cell 

36 

62.08 

2.53 

To  validate  the  FE  modeling  technique  further,  data  on  plain  weave  composites  available  in 
literature  is  considered.  Three  balanced  plain  weaves,  viz.  E-Glass/Epoxy  [10], 
AS4-Carbon/Epoxy  [9],  and  Carbon/Epoxy  plain  woven  fabrics  [6],  are  studied.  The  boundary 
conditions  described  in  Eqn.  1  and  Eqn.  2  are  used  to  determine  the  elastic  moduli  of  plain 
weaves.  Table  5  lists  the  material  properties  of  the  constituents  of  the  weaves.  The  geometri¬ 
cal  properties  of  the  woven  fabrics  used  in  the  computations  are  shown  in  Table  6.  A  compari¬ 
son  of  the  elastic  constants  of  plain  weaves  obtained  from  the  present  FE  model  with  analytical 
and  experimental  results  is  shown  in  Table  7.  As  can  be  observed,  there  is  a  good  correlation 
between  the  present  resuits  and  those  available  in  literature.  The  slight  discrepancies  in  vaiues 
arise  as  a  result  of  differences  in  modeling  techniques  between  the  present  approach  and 
those  adopted  in  literature. 


Table  5.  Elastic  Properties  of  Yams  and  Matrix  Used  in  the  Analysis 


Material 

Material  Property 

Ell 

E22 

Vl2 

Gi2 

G23 

■bh 

■Bl 

MBM 

mmim 

E-Glass/Epoxy  [10] 
(V/=0.65) 

48.47 

(7.03) 

18.06 

(2.62) 

0.25 

0.34 

5.58 

(0.81) 

3.31 

(0.48) 

AS4-Carbon/Epoxy  [9] 
(V/=  0.70) 

155.83 

(22.60) 

10.13 

(1.47) 

0.24 

0.5 

5.72 

(0.83) 

3.38 

(0.49) 

Carbon/Epoxy  [6] 

(¥/=  0.6) 

135.27 

(19.61) 

9.65 

(1.40) 

0.28 

0.43 

5.37 

(0.78) 

3.38 

(0.49) 

Epoxy  3502  [9] 

mswEm 

mSmm 

0.35 

0.35 

■m 

mBm 
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Table  6.  Geometrical  Properties  of  Yams 


Material 

MSSM 

Yarn  thickness 
mm  (in.) 

Yarn  spacing 
mm  (in.) 

V/ 

Vf" 

E-Glass/epoxy  [10] 

1.747 

(0.068) 

0.268 

(0.010) 

1.747 

(0.068) 

0.65 

0.35 

AS4-Carbon/epoxy  [9] 

1.760 

(0.069) 

0.129 

(0.005) 

1.760 

(0.069) 

0.70 

0.60 

Carbon/epoxy  [6] 

2.499 

(0.0984) 

0.109 

(0.0043) 

2.499 

(0.0984) 

0.63 

0.50 

Table  7.  Comparison  of  Present  FEA  Results  With  Other  Analytical/Experimentai 
Results  for  Plain  Woven  Fabrics 


Material 

Property 

Present 

KBH 

E-Glass/epoxy 
plain  weave  [1 0] 
Vf®=035 

E,  GPa  (Msi) 

Vxy 

0.14 

IHIHI 

Gxy  GPa  (Msi) 

5.972 

(0.866) 

* 

AS4-Carbon/epoxy 
plain  weave  [9] 

Vf  =0.60 

Ex  GPa  (Msi) 

60.40 

(8.760) 

60.0 

(8.702) 

^xy 

0.043 

- 

II^HH 

Gxy  GPa  (Msi) 

■jjEKglH 

Carbon/epoxy  plain 
weave  [6] 

Vf“=030 

Ex  GPa  (Msi) 

26.1 

(3.785) 

25.0 

(3.62) 

Vxy 

0.06 

- 

Gxy  GPa  (Msi) 

wmassam 

The  reduction  in  the  in-piane  stiffness  of  woven  fabrics  in  comparison  with  that  of  unidirec¬ 
tional  lamina  is  infiuenced  to  a  large  extent  by  the  waviness  of  the  yarns.  This  degree  of  un¬ 
dulation  of  the  yarns  depends  on  the  thickness  (h)  and  width  (a)  of  the  yarns.  Lower  h/a  ratio 
indicates  a  lower  degree  of  undulation  and  vice  versa.  The  yarn  fiber  volume  fraction  is  kept 
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constant  (V/  =  0.70).  It  can  be  seen  that  the  modulus  E,  reduces  with  increase  in  h/a  values. 
As  h/a  values  increases,  so  does  unduiation.  This  is  accompanied  by  a  iowering  of  the  overaii 
volume  fraction  (Vf°).  However,  the  Poisson’s  ratio  decreases  with  increase  in  overall  volume 
fraction.  Table  8  lists  the  variation  of  elastic  properties  with  the  degree  of  undulation  (h/a)  and 
the  overall  fiber  volume  fraction  of  the  weave.  It  can  be  observed  that  as  the  undulation  de¬ 
creases,  the  yarns  get  compacted  leading  to  an  increase  in  the  overall  fiber  volume  fraction. 
With  increase  in  undulation  length,  the  magnitude  of  the  local  fiber  angle  decreases.  This  also 
reduces  the  coupling  between  extension  and  bending  effects.  The  shear  modulus  increases 
with  increase  in  the  degree  of  undulation  (h/a)  of  the  yarns. 


Table  8.  Effect  of  h/a  and  u/a  on  Material  Properties  of  AS4-Carbon/Epoxy 
(V/  =  0.70)  Plain  Weave  [9] 


Vf" 

h/a 

■Ml 

0.45 

0.172 

(0.006) 

1.76 

(0.069) 

0.098 

46.13 

(6.69) 

0.053 

0.50 

0.154 

(0.006) 

1.76 

(0.069) 

0.088 

0.049 

IHI 

0.55 

0.143 

(0.005) 

1.76 

(0.069) 

0.082 

mSmm 

0.046 

■EBBI 

0.60 

0.129 

(0.005) 

1.76 

(0.069) 

0.074 

60.40 

(8.76) 

0.044 

4.883 

(0.708) 

0.65 

Mi 

0.068 

0.041 

4.854 

(0.703) 

0.70 

0.111 

(0.004) 

1.76 

(0.069) 

0.064 

72.40 

(10-5) 

0.036 

4.863 

(0.705) 

FABRIC  STRESS  DISTRIBUTIONS 


A  study  of  the  internal  stress  distributions  generated  within  the  fabric  when  subjected  to 
external  loads  is  carried  out.  This  wouid  be  useful  in  identifying  the  dominant  ones  and  also 
regions  of  stress  concentrations.  This  information  would  aid  in  predicting  possibie  modes  of 
damage  likely  to  occur  within  the  yarns  and  surrounding  pure  matrix  pockets.  Such  a  study 
wouid  also  help  in  assessing  the  applicabiiity  of  damage  models  for  describing  the  nonlinear 
behavior  of  yarn  and  pure  matrix  pockets. 
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For  the  sake  of  conciseness,  only  the  important  stresses  generated  within  piain  weave  fab¬ 
rics  as  a  result  of  tensile  loading  are  discussed  here.  The  carbon/epoxy  piain  weave  (Vf®  =  0.5) 
considered  eariier  is  subjected  to  tension  aiong  the  warp  direction.  Riots  of  the  stresses  gener¬ 
ated  in  the  warp  and  fill  corresponding  to  an  overall  applied  stress  an  =  1 83.7  MPa  (26.64  ksi) 
are  shown  in  Figures  6  - 1 1 .  All  the  yarn  stresses  are  measured  with  respect  to  the  local  mate¬ 
rial  axes.  It  can  be  seen  that  at  regions  close  to  the  point  where  there  is  a  change  in  curvature 
of  the  warp,  high  values  of  in-plane  shear  stress  and  on  are  prevalent.  Such  stress  concen¬ 
trations  can  trigger  localized  damage  in  the  warp  even  when  the  overall  applied  stress  (an)  is 
much  below  the  final  strength  of  the  weave.  The  out-of-plane  normal  stresses  033  developed  in 
the  warp,  though  smaller  in  comparison  with  the  in-plane  stresses  an  (peak  033  =  4.682  MPa 
(0.679  ksi)  against  a  peak  aii=  454.64  MPa  (65.938  ksi))  can  cause  some  damage  in  the 
yarns.  The  variation  of  the  in-plane  normal  stress  (an)  across  the  thickness  of  the  warp  indi¬ 
cates  localized  bending  occurring  as  a  result  of  extension.  Also,  small  regions  can  be  seen 
where  the  local  in-plane  stresses  (an  and  x^)  are  very  high.  The  transverse  normal  stress  (022) 
in  the  fill  (orthogonal  to  the  direction  of  loading)  is  higher  than  that  in  the  warp  and  hence 
transverse  tensile  failure  in  fill  would  occur  much  before  final  failure  of  the  weave. 


Figure  6.  Tensile  Stress  (an)  Developed  in  Warp  of  Carbon/Epoxy  (Vf  =  0.5) 
Subjected  to  Tension  Along  Warp  Corresponding  to  Overall  Applied 
Tensile  Stress  of  183.7  MPa. 

(All  stresses  in  MPa  and  are  measured  with  reference  to  the  local 
material  axes) 
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Figure  7.  Out-of-Plane  Stress  (033)  Developed  in  Warp  of  Carbon/Epoxy  (Vf 
0.5)  Subjected  to  Tension  Along  Warp  Corresponding  to  Applied 
Overall  Tensile  Stress  of  183.7MPa. 

(All  stresses  in  MPa  and  are  measured  with  reference  to  the  local 
material  axes) 


Figure  8.  In-Plane  Shear  Stress  (012)  Developed  in  warp  of  Carbon/Epoxy  (Vf 
0.5)  Subjected  to  Overall  Tensile  Stress  of  183.65  MPa  (26.64  ksi). 
(All  stresses  in  MPa  and  are  measured  with  reference  to  the  local 
material  axes) 


Figure  10.  Transverse  Normal  Stress  (022)  Developed  in  Fill  of  Carbon/Epoxy  (Vf 
=  0.5)  Subjected  to  Overall  Tensile  Stress  of  183.65  MPa. 

(All  stresses  in  MPa  and  are  measured  with  reference  to  the  local 
material  axes) 
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Figure  11.  Transverse  Shear  Stress  ('T23)  Developed  in  Fill  of  Carbon/Epoxy  (Vf= 
0.5)  Subjected  to  Overall  Tensile  Stress  of  183.65  MPa  (26.64  ksi). 
(All  stresses  in  MPa  and  are  measured  with  reference  to  the  local 
material  axes) 

MATERIAL  MODELING 


The  nonlinear  behavior  of  the  weave  constituents,  viz.  yarns  and  pure  matrix,  is  modeled 
using  a  damage  mechanics  based  approach.  The  dissipated  energy  density  (<]))  is  chosen  as 
the  damage  parameter.  As  discussed  in  Appendix  B,  the  nature  of  a  damage  surface  for  a 
fixed  value  of  ^  and  its  change  with  ^  can  be  modeled  at  various  levels  of  complexity  to  simu¬ 
late  the  “in-situ”  nonlinear  material  behaviors.  In  the  results  reported  here,  simple  quadratic 
forms  of  the  damage  surfaces  in  terms  of  chosen  strain  variables  are  employed.  The  changes 
of  these  surfaces  with  (j)  are  chosen  in  a  manner  such  that  they  yield  power  law  type  stress- 
strain  relations  under  unidirectional  straining  beyond  the  initial  elastic  domain  as  described  in 
Appendix  B.  The  material  characterization  curves  for  AS4-Carbon/Epoxy  and  E-Glass/Epoxy 
unidirectional  lamina  predicted  by  the  current  methodology  compare  well  with  experimental  re¬ 
sults  contained  in  References  [15]  and  [16],  respectively.  The  nonlinear  behavior  of  the  textile 
constituents  used  in  the  present  analyses  is  shown  in  Figures  12-17. 

The  material  constitutive  laws  for  yams  and  pure  matrix  described  in  Appendix  B  are  incor¬ 
porated  as  a  user  subroutine  in  ABACUS  [13].  The  material  properties  required  for  analyses 
are  the  elastic  moduli,  Poisson’s  ratios  and  variables  that  describe  the  damage  surfaces  in  the 
strain  space  for  various  values  of  dissipated  energy  density  <|)  in  the  yarn/matrix. 
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Figure  12.  Longitudinal  Stress-Strain  Curve  for  Epoxy  3502 


Shear  strain 


Figure  1 3.  Shear  Stress-Strain  Curve  for  Epoxy  3502 
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Shear  stress  (MPa) 


Figure  14.  Longitudinal  Stress-Strain  Curve  for  AS4/Epoxy  (Vf  =  0.7) 


Shear  strain 


Figure  15.  Shear  Stress-Strain  Curve  for  AS4/Epoxy  (Vf  =  0.7) 
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NONLINEAR  ANALYSIS 


A  typical  unit  cell  of  plain  weave  discretized  by  finite  elements  is  shown  in  Figure  3.  Within 
each  finite  element  of  the  unit  cell,  the  stiffnesses,  stresses  and  strains  can  be  monitored  at 
the  material  integration  points.  This  function  is  performed  by  the  material  models  described  in 
Appendix  B.  For  this  purpose,  the  constitutive  laws  for  the  materials  are  implemented  through  a 
user-defined  subroutine  (UMAT)  and  linked  with  the  finite  element  analysis  program,  ABAQUS 
[13].  During  each  load  increment  of  the  nonlinear  analysis  process,  the  material  definition  rou¬ 
tines  describe  the  current  stiffnesses  and  damage  state  at  the  numerical  integration  points  of 
the  finite  elements.  The  reductions  in  stiffnesses  occurring  as  a  result  of  localized  damage 
within  the  yarn/pure  matrix  pockets  are  monitored  and  incorporated  into  the  weave  response. 

An  incremental-iterative  approach  is  adopted  for  the  nonlinear  finite  element  analysis.  The 
modified  Newton-Raphson  method  is  used  to  trace  the  loading  path  of  the  structure.  The  non¬ 
linear  analysis  starts  from  an  equilibrium  state.  Loads  are  applied  by  prescribing  incremental 
displacements  and  the  unbalanced  forces  evaluated.  An  iterative  procedure  is  used  until  the 
out-of-balance  load  vector  or  the  displacement  increments  are  sufficiently  small.  Within  each 
element  representing  yarn/pure  matrix,  each  Gaussian  integration  point  represents  a  certain 
volume  of  material  whose  material  stiffnesses  are  described  by  the  user  material  definition 
subroutine  developed.  Thus,  for  a  particular  increment  of  load,  depending  on  the  current  value 
of  the  dissipated  energy  density  (|),  an  element  may  contain  some  points  with  no  damage  and 
some  with  damage.  Within  each  load  increment  of  the  nonlinear  analysis,  the  material  is 
checked  for  occurrence  of  damage.  Once  damage  is  detected,  the  material  properties  are 
changed  accordingly.  At  any  given  stage  of  loading,  the  above  approach  of  detecting  damage 
is  applied  at  all  Gaussian  integration  points  of  the  yarn/pure  matrix  elements.  Thus,  by  model¬ 
ing  damage  at  discrete  integration  points,  the  progressive  failure  of  the  woven  fabric  is  ana¬ 
lyzed. 

The  behavior  of  balanced  E-Glass/Epoxy  [10]  and  AS4-Carbon/Epoxy  [9]  and  car¬ 
bon/epoxy  [6]  plain  woven  fabric,  described  earlier,  when  subjected  to  in-plane  loads  is  exam¬ 
ined.  The  current  predictions  are  compared  with  experimental  results  if  available  or  with  other 
theoretical  results  available  in  literature.  Elastic  properties  of  the  weave  constituents  are  tabu¬ 
lated  in  Table  5.  The  geometrical  parameters  of  the  weave  used  in  the  present  investigation  are 
listed  in  Table  6. 
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TENSILE  BEHAVIOR 


The  plain  woven  fabrics  are  subjected  to  tensile  loading  along  warp  using  boundary  condi¬ 
tions  described  in  Eqn.  1.  Effects  of  geometric  nonlinearity  and  materiai  nonlinearity  of  the 
yarns  and  surrounding  pure  matrix  are  included  in  the  analysis.  Figures  18  and  19  show  the 
stress/strain  curves  obtained  from  the  present  FE  analysis. 

These  results  are  compared  with  experimentai  values  [6,  9,  10].  Figure  18  shows  a  good 
agreement  between  the  resuits  of  the  present  analysis  and  experimentai  values  in  the  pre- 
’knee’  portion  of  the  stress-strain  curve.  The  ‘knee’  in  the  stress-strain  curve  is  caused  due  to 
transverse  tensile  failure  in  the  fill. 

Carbon/epoxy  unidirectional  lamina  is  stiffer  than  E-Glass/Epoxy  lamina.  As  a  result,  inter¬ 
mediate  damage  mechanisms  (e.g.,  ‘knee’)  prior  to  fracture  of  fiber  are  not  seen  in  the 
stress/strain  response  of  AS4-Carbon/Epoxy  plain  weave  (Figure  19)  as  clearly  as  that  in 
E-Glass/epoxy  piain  weave  (Figure  18).  Hence,  the  stress/strain  plot  is  almost  linear  until  final 
failure.  On  the  other  hand,  in  the  case  of  E-Glass/Epoxy  plain  weave,  the  ‘knee’  is  more  pro¬ 
nounced  than  that  in  AS4-Carbon/Epoxy  plain  weave.  The  damage  mechanisms  seen  in  the 
plain  weaves  subjected  to  tension  aiong  warp  direction  are;  (i)  Transverse  tensile  failure  (of 
matrix)  in  fiii;  (ii)  Matrix  failure  in  in-piane  shear  in  fiii  near  the  yarn  extremities  where  there  is  a 
stress  concentration;  and  (iii)  Fiber  failure  in  tension  in  warp. 

Figure  20  shows  the  behavior  of  carbon/epoxy  plain  weave  lamina  [6]  subjected  to  tension 
aiong  warp  and  tension  along  an  off-axis  (0  =  IS”)  direction.  The  progression  of  damage  within 
the  fabric  when  the  fabric  is  subjected  to  tension  along  warp  is  shown  by  contour  piots  of  the 
damage  parameter  ^  in  Figures  21-25.  For  the  case  of  off-axis  loading,  considerable  shear 
softening  occurs  in  the  transverse  fill  before  final  failure  of  the  weave.  Damage  in  the  fiii  due  to 
transverse  tension  occurs  much  earlier  (177  MPa)  than  the  damage  in  the  warp  caused  by  fi¬ 
ber  failure.  The  damage  mechanisms  seen  in  the  plain  weave  subjected  to  off-axis  tension  are: 
(i)  Transverse  tensile  failure  (of  matrix)  in  fiii;  (ii)  Matrix  faiiure  in  shear  in  fill;  and  (iii)  Fiber  faiiure 
in  tension  in  warp. 
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Figure  18.  Response  of  E-Glass/Epoxy  Plain  Weave  Subjected  to  Tension 


Figure  19.  Response  of  AS4/Carbon/Epoxy  Plain  Weave  Subjected  to  Tension 
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A  Experimental  [8]  (15  degree) 
X  Experimental  [8]  (0  degree) 

- Present  (0  degree) 

Present  (15  degrees) _ 


Twistle  Strain 


Figure  20:  Response  of  Carbon/Epoxy  [6]  Plain  Weave  Subjected  to  Tension 


Figure  21 .  Variation  of  ^  in  Fill  of  Carbon/Epoxy  Plain  Weave  Subjected  to 
Tension  Aiong  Warp  (Oo  =168.10  Mpa ) 


Figure  22.  Variation  of  <l)  in  Fill  of  Carbon/Epoxy  Plain  Weave  Subjected  to 
Tension  Along  Warp  (Oo  =  188.9  MPa) 


Figure  23.  Variation  of  <|)  in  Fill  of  Carbon/Epoxy  Plain  Weave  Subjected  to 
Tension  Along  Warp  (ao=  210.05  MPa) 


Figure  24.  Variation  of  ^  in  Warp  of  Carbon/Epoxy  Plain  Weave  Subjected  to 
Tension  Along  Warp  (Oo  =  421 .63  MPa ) 
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Figure  25.  Variation  of  ^  in  Warp  of  Carbon/Epoxy  Piain  Weave  Subjected  to 
Tension  Along  Warp  (Oo  =  433.12  MPa) 


IN-PLANE  SHEAR 

The  behavior  of  E-Glass/Epoxy  and  AS4-Carbon/Epoxy  plain  weaves  under  in-plane  shear 
is  examined  using  boundary  conditions  shown  in  Eqn.  2.  Figures  20  and  21  show  the  stress- 
strain  curves  obtained  from  the  present  analysis.  For  the  case  of  AS4-Carbon/Epoxy  plain 
weave,  the  present  results  are  compared  with  an  experimental  result  available  in  literature  [9]. 
The  ultimate  strength  prediction  agrees  well  with  the  experimental  result  though  the  failure 
strains  differ.  This  could  be  due  to  differences  in  material  properties  used  in  the  present  mate¬ 
rial  model  and  those  tested.  The  sequence  of  failure  for  plain  weaves  subjected  to  in-plane 
shear  is:  (i)  Matrix  failure  in  shear  in  fill;  (ii)  Matrix  failure  in  out-of-plane  shear  in  warp;  and  (iii) 
Matrix  failure  in  transverse  tension  in  fill. 


Strms(MPa) 


CONCLUSIONS 


An  investigation  into  the  micromechanics  of  fabric  composites  requires  a  proper  description 
of  the  fabric  architecture  and  varying  material  properties.  This  has  been  accomplished  by  the 
use  of  3D  finite  elements  which  permits  a  detailed  modeling  of  the  varying  cross-sectional  pro¬ 
file  and  undulations  of  the  yarns.  The  mesh  generation  routines  which  have  been  developed  in 
the  current  study  can  be  used  to  generate  FE  models  for  a  variety  of  woven  fabrics.  The  output 
from  the  mesh  generation  programs  serve  as  input  files  for  a  general-purpose  finite  element 
analysis  program,  ABACUS. 

The  periodic  nature  of  the  weave  pattern  combined  with  point  symmetry  of  material  distri¬ 
bution  can  be  used  advantageously  in  the  analysis.  This  results  in  a  reduction  of  the  problem 
size  and  enhances  computational  efficiency. 

A  study  of  the  stress  distributions  generated  within  the  fabrics  indicates  the  presence  of 
concentrations  of  high  values  of  stress  in  regions  with  sharp  changes  in  yarn  curvature.  This 
would  lead  to  localized  damage  even  at  stresses  much  below  the  ultimate  strength  of  the 
weave.  Hence,  material  models  which  monitor  the  onset  of  damage  and  accordingly  reflect 
them  appropriately  in  the  properties  of  the  fabric  constituents  are  necessary. 

Damage  mechanics  based  constitutive  laws  based  on  strain  energy  dissipation  (|>  as  the 
damage  parameter  is  developed.  The  nature  of  damage  surface  for  a  particular  value  of  (j)  and 
its  change  with  ^  is  modeled  to  simulate  the  “in-situ”  nonlinear  material  behaviors.  In  the  pres¬ 
ent  study,  simple  quadratic  forms  of  the  damage  surfaces  in  terms  of  chosen  strain  variables 
are  employed.  The  changes  of  these  surfaces  with  ^  are  chosen  in  a  manner  such  that  they 
yield  power  law  type  stress-strain  relations  under  unidirectional  straining  beyond  the  initial 
elastic  domain.  The  material  model  has  been  used  to  study  and  describe  the  initiation  and 
growth  of  damage  in  plain  woven  fabrics.  The  damage  modeling  methodology  presented  herein 
has  wide  applicability  in  describing  response  of  materials  which  do  not  show  any  rate  depend¬ 
ent  or  viscous  behavior. 

A  progressive  failure  analysis  of  plain  woven  fabrics  subjected  to  in-plane  loading  has  been 
undertaken  and  the  results  compared  with  analytical/experimental  results  in  literature.  Both 
material  and  geometric  nonlinearities  are  considered.  Plain  weaves,  when  subjected  to  tension, 
exhibit  an  almost  linear  behavior  prior  to  failure.  The  dominant  failure  mechanism  is  transverse 
tensile  failure  of  the  fill  for  a  plain  weave  loaded  in  the  warp  direction.  Such  failure  causes  a 
reduction  in  the  stiffness  of  the  plain  weave  and  is  reflected  as  a  ‘knee’  in  the  stress-strain  plot. 
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Plain  woven  fabrics  exhibit  a  nonlinear  behavior  when  subjected  to  in-plane  shear.  This  nonlin¬ 
ear  behavior  is  mainly  due  to  the  shear  softening  of  matrix  in  the  yarns.  Failure  of  matrix  in 
warp  due  to  in-plane  shear  and  transverse  tensile  failure  of  the  fill  are  prime  mechanisms  of 
failure  under  such  loading  conditions. 
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INTERACTIONS/TRANSITIONS 
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APPENDIX  A  -  BOUNDARY  CONDITIONS 


The  boundary  conditions  for  a  periodic  cell  of  balanced  woven  fabric  are  first  stated.  The 
woven  fabric  is  assumed  to  be  a  part  of  a  laminate  in  which  successive  woven  fabrics  are 
stacked  symmetricaiiy  about  the  xy  plane. 


Figure  A1 .  Periodic  Cell 


The  conditions  which  need  to  be  satisfied  on  account  of  the  periodic  nature  of  the  weave  ge¬ 
ometry  are  (Figure  A1 ): 


On  the  planes  x  =  ±  b, 

u(b,y,z)  -  u(-b,y,z)  =  ua  -  uc 
v(b,y,z)  -  v(-b,y,z)  =  va  -  vc 
w(b,y,z)  -  w(-b,y,z)  =  wa  -  wc 

On  the  planes  y  =  ±  b, 

U(x,b,z)  -  U(x,-b,z)  =  UA  -  UB 
V(x,b,z)  -  V(x,-b,z)  =  VA  -  VB 
w(x,b,z)  -  w(x,-b,z)  =  WA  -  WB 


(A1) 


(A2) 
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(A3) 


On  the  planes  z  =  ±t, 

u(x,y,t)  -  u(x,y,-t)  =  ua  *  ud 
v(x,y,t)  -  v(x,y,-t)  =  va  -  vd 
w(x,y,t)  -  w(x,y,-t)  =  wa  -  wd 


A1.1  BOUNDARY  CONDITIONS  FOR  IN-PLANE  EXTENSION 


A1.1.1  Boundary  Conditions  for  Unit  Cell 


An  extension  along  the  x-axis  is  prescribed  such  that 


UA  = 

-UC 

=  Uo 

VA  = 

-VB 

==  Vo 

WA  = 

“  WD 

=  Wo 

The  periodic  conditions  for  extension  from  Eqns  A1  -A3  are: 

On  the  planes  x  =  ±b, 

u(b,y,z)-u(-b,y,z)  =  2uo 

v(b,y,z)  -  v(-b,y,z)  =  0  (A4) 

w(b,y,z)  -  w(-b,y,z)  =  0 

On  the  planes  y  =  ±b, 

u(x,b,z)  -  u(x,-b,z)  =  0 

v(x,b,z)  -  v(x,-b,z)  =  2  Vo  (A5) 

w(x,b,z)  -  w(x,-b,z)  =  0 

On  the  planes  w  =  ±t, 

u(x,y,t)  -  u(x,y,-t)  =  0 

v(x,y,t)  -  v(x,y,-t)  =  0  (A6) 

w(x,y,t)  -  w(x,y,-t)  =  2wo 
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Figure  A2.  Unit  Cell  of  Plain  Weave 

Now,  for  symmetry  about  the  xy,  yz,  xz  planes,  the  following  conditions  hold  good: 

About  xy  plane: 

u(x,y,z)  =  u(x,y,-z) 

v(x,y,z)  =  v(x,y,-z)  (A7) 

w(x,y,z)  =  -w(x,y,-z) 

About  the  xz  plane: 

u(x,y,z)  =  u(x,-y,z) 

v(x,y,z)  =  -v(x,-y,z)  (A8) 

w(x,y,z)  =  w(x,-y,z) 

About  the  yz  plane: 

u(x,y,z)  =  -u(-x,y,z) 

v(x,y,z)  =  v(-x,y,z)  (A9) 

w(x,y,z)  =  w(-x,y,z) 

Applying  these  (Eqn.  A7-A9)  symmetries  to  Eqn.  A4  -  A6,  the  boundary  conditions  (Figure  A2) 
are: 
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(A10) 


u(-,y,z)  =  -u(--,y,z)  =  uo 

v(x,|,z)  =  -v(x,-Y>z)  =  Vo 
w(x,y,t)  =  -w(x,y,-t)  =  Wo 


A1.1.2  Boundary  Conditions  for  Quarter  Cell 

The  antisymmetry  of  material  distribution  about  the  planes  x  =  0  and  symmetry  about  the 
planes  x  =  ±  b/2  can  be  used  to  derive  conditions  applicable  to  each  quarter  of  the  unit  cell  (Fig¬ 
ure  A3). 


Figure  A3.  Quarter  of  Plain  Weave  Unit  Cell 

For  extension  in  the  x  -  direction,  from  the  antisymmetry  present  in  the  weave,  we  note  that 
on  a  x-plane,  the  following  conditions  are  valid: 

u(x,y,z)=  u(x,-y,-z) 

about  the  X  -  axis  :  v(x,y,z)  =  -v(x,-y,-z)  (A1 1 ) 

w(x,y,z)  =  -w(x,-y,-z) 
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Similarly,  the  anti-symmetric  stiffness  distribution  about  the  y-axis  gives  rise  to  the  following 
conditions: 


u(x,y,z)  =  -u(-x,y,-z) 

about  the  Y  axis  :  v(x,y,z)=  v(-x,y,-z)  (A12) 

w(x,y,z)  =  -w(-x,y,-z) 


Finally,  there  is  point-symmetry  about  the  z-axis.  This  gives  us  the  relations: 

u(x,y,z)  =  -u(-x,-y,z) 

about  the  Z  axis :  v(x,y,z)  =  -v(-x,-y,z)  (A1 3) 

w(x,y,z)  =  w(-x,-y,z) 

Eqns.  A10-A12  can  now  be  used  to  obtain  the  relationship  between  points  lying  on  the  xy,  zx 
and  yz  planes. 


u(x,y,0)  =  -u(-x,y,0)  =  u(x,-y,0) 
on  the  xy  plane(z  =  0)  :  v(x,y,0)  =  v(-x,y,0)  =  -v(x,-y,0) 

w(x,y,0)  =  w(-x,y,0)  =  -w(x,-y,0) 
u(0,y,z)  =  -u(0,y,-z)  =  -u(0,-y,z) 

on  the  yz  plane(x  =  0) :  v(0,y,z)  =  v(0,y,-z)  =  -v(0,-y,z)  (A1 4) 

w(0,y,z)  =  -w(0,y,-z)  =  w(0,-y,z) 
u(x,0,z)  =  -u(-x,0,z)  =  u(x,0,-z) 
on  the  xz  plane  (y  =  0) :  v(x,0,z)  =  v(-x,0,z)  =  -v(x,0,z) 

w(x,0,z)  =  w(-x,0,z)  =  -w(x,0,-z) 

Eqn.  A10  specifies  the  boundary  conditions  for  an  unit  ce// with  respect  to  the  coordinate  axes 
shown  in  Figure  1 .  With  respect  to  a  coordinate  axes  having  its  origin  at  the  center  of  the  unit 
cell  (in  the  X  >  0,  y  >  0  quadrant),  we  get  the  boundary  conditions  for  a  unit  cell  subjected  to  ex¬ 
tension  as  (Figure  2): 
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Ontheplanex  =  0  : 

u(0,y,z)  =  -u(0,y,-z) 
v(0,y,z)  =  v(0,y,-z) 
w(0,y,z)  =  -w(0,y,-z) 

Alongthey  -  axis  (0,y,0)  ; 

u(0,y,0)  =  0.0 

Ontheplaney  =  0  ; 

u(x,0,z)  =  u(x,0,-z) 
v(x,0,z)  =  -v(x,0,-z) 
w(x,0,z)  =  -w(x,0,-z) 

Alongthex-axis  (x,0,0) 

;  v(x,0,0)  =  0.0 

Along  thez-  axis  (0,0,  z) : 

u(0,0,z)  =  0 
v(0,0,z)  =  0 
w(0,0,z)  =  -w(0,0,-z) 

On  the  plane  x  =  b/2  ; 

.b  uo 

"<!’>'•==)  y 

Ontheplaney=b/2: 

v(x,y,z)  =  uniform 

A1.2  BOUNDARY  CONDITIONS  FOR  IN-PLANE  SHEAR 


The  periodic  cell  is  subjected  to  in-plane  shear  loading  using  the  following  conditions: 


UA  -  -UB  “  P 
VA  =  -VC  =  P 
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Substituting  the  above  into  Eqn.  (A1-A3),  we  get  the  following  periodic  conditions  for  in¬ 
plane  shear; 

Along  the  planes  x  =  ±b, 

u(b,y,z)  -  u(-b,y,z)  =  0 

v(b,y,z)  -  v(-b,y,z)  =  2p  (A1 7) 

w(b,y,z)-w(-b,y,z)  =  0 

Along  the  planes  y  =  ±b, 

u(x,b,z)  -  u(x,-b,z)  =  2p 

v(x,b,z)  -  v(x,-b,z)  =  0  (A1 8) 

w(x,b,z)  -  w(x,-b,z)  =  0 

Along  the  planes  z  =  ±  t, 

u(x,y,t)  -  u(x,y,-t)  =  0 

v(x,y,t)  -  v(x,y,-t)  =  0  (A1 9) 

w(x,y,t)  -  w(x,y,-t)  =  0 

Now,  about  the  xy  plane,  we  have, 
u(x,y,z)  =  u(x,y,-z) 

v(x,y,z)  =  v(x,y,-z)  (A20) 

w(x,y,z)  =  -w(x,y,-z) 

About  the  xz  plane, 

u(x,y,z)  =  -u(x,-y,z) 

v(x,y,z)  =  v(x,-y,z)  (A21) 

w(x,y,z)  =  -w(x,-y,z) 

About  the  yz  plane, 

u(x,y,z)  =  u(-x,y,z) 

v(x,y,z)  =  -v(-x,y,z)  (A22) 

w(x,y,z)  =  -w(-x,y,z) 

Using  Eqn.  A20,  we  get 
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w(x,y,0)  =  0 


Also  at  z  =  ±t,  from  Eqn  A1 9  and  Eqn.  A20  we  have 
w(x,y,t)  =  w(x,y,-t)  =  0 


From  Eqn.  A22,  at  x  =  0,  we  get 


v(0,y,z)  =  0 
w(0,y,z)  =  0 


From  Eqn  A21 ,  at  y  =  0,  we  get 


u(x,0,z)  =  0 
w(x,0,z)  =  0 

Eqn.  A20-A22  together  with  the  periodic  conditions,  Eqn.  A17-A19,  yieid  the  foiiowing  condi¬ 
tions: 

Eqn  A17&A20  =>  u(x,b,z)  =  -u(x,-b,z)  =  p 
Eqn.  A16&A21  =>  v(b,y,z)  =  -v(-b,y,z)  =  p 

Eqn.A16&A21  =>  w(±b,y,z)  =  0  (A23) 

Eqn.  A 1 7  &  A20  =>  w(x,  ±  b,z)  =  0 
Eqn.  A1 8&  A 1 9  w(x,y,  + 1)  =  0 

A1.2.2  Boundary  Conditions  for  Quarter  Cell 

Boundary  conditions  for  a  quarter  celi  (Figure  3)  can  now  be  stated  using  those  of  the  unit 
cell.  Antisymmetries  of  loading  about  the  x  and  y  axis  give  the  foiiowing  reiations; 


A-8 


About  the  x  -  axis, 


u(x,y,z)  =  -u(x,-y,-z) 

v(x,y,z)  =  v(x,-y,-z)  (A24) 

w(x,y,z)  =  w(x,-y,-z) 

About  the  y  -  axis, 

u(x,y,z)=  u(-x,y,-z) 

v(x,y,z)  =  -v(-x,y,-z)  (A25) 

w(x,y,z)  =  w(-x,y,-z) 

About  the  z  -  axis 

u(x,y,z)  =  -u(-x,-y,z) 

v(x,y,z)  =  -v(-x,-y,z)  (A26) 

w(x,y,z)  =  w(-x,-y,z) 

Eqn.  (A23)  and  Eqns.  (A24-A26)  can  be  used  to  obtain  the  boundary  conditions  prevalent  for 
the  quarter  of  unit  cell  shown  in  Figure  A3  as: 

On  the  plane  x  =  0  u(0,y,z)  =  u(0,y,-z) 
v(0,y,z)  =  -v(0,y,-z) 
w(0,y,z)  =  w(0,y,-z) 

On  the  plane  y  =  0  u(x,0,z)  =  -u(x,0,-z) 
v(x,0,z)  =  v(x,0,-z) 
w(x,0,z)  =  w(x,0,-z) 

Ontheplaney=— 


Ontheplanex  =  — 


u(x,|,z)=£ 
w(x,-,z)  =  0 


w(-,y,z)  =  0 


(A27) 
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APPENDIX  B  -  DAMAGE  MECHANICS  BASED  CONSTITUTIVE  LAWS 


BACKGROUND 


Constitutive  laws  for  continuum  damage  mechanics  are  discussed  in  [B1-B3].  Such  models 
have  been  found  to  be  very  successful  representing  the  behavior  of  brittle  materials  and  make 
use  of  a  phenomenological  approach  and  the  internal  variables  o",  which  are  the  damage  vari¬ 
ables  (such  as  average  crack  densities).  One  can  write  the  formulation  given  in  Reference  2  in 
the  following  form  where  a  and  s  are  vectors. 


g  •  ds  -  pdvj/  =d(j)  >  0;  Clausius-Duhem  inequality  (see  Figure  A1 )  (B1  a) 

p  =  mass  density 

g  =  g[8,  (B“;a=  1,2,  ...k] 

k  =  total  number  of  damage  variables 

v|/  =  Helmholtz  free  energy 

=  Vj/  [8,  <0“] 

(B1b) 

^£i 


(|)  =  <|)  (8,  co“),  the  dissipated  energy  density 

d<!>=  i  R“dfi)®^0  (B1c) 

0=1 

R“  =  conjugate  thermodynamic  force,  an  energy  release  rate  for  the 
damage  variable  ©“  =  R“  (s,  (o“) 


3/r  _  d<l> 


(Bid) 


In  the  strain  energy  dissipation  (SED)  concept  employed  in  [B4],  there  is  one  damage  vari¬ 
able  (k=  1)  which  is  the  dissipated  energy  density  (|)  (=(o’),  and  for  this  case  one  may  choose  the 


following  form  for  the  Helmholtz  free  energy 
PW = 

^  O 

where  <j)c  is  the  current  value  of  (j)  and  Qj  are  the  elastic  stiffnesses  in  contracted  notation.  It  fol¬ 
lows  from  the  expression  for  R*  =  R  in  (Bid)  that  for  damage  growth  to  occur 

R(f,^)  =  l  (B3) 
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Figure  A1 .  Helmholtz  Free  Energy  pv|/  and  Dissipated  Energy  (]) 

which  defines  the  damage  surface  (for  a  given  value  of  (|))  in  the  strain  space  and  when 
R  <  1 ,  no  damage  growth  would  occur.  It  is  clear  that  the  dissipated  energy  is  a  function 
of  the  strain  variables  independent  of  the  loading  path  and  unloading  will  occur  in  a  straight  line 
from  the  current  strain  state  to  the  origin.  Obviously,  these  assumptions  have  to  be  checked  by 
conducting  tests  and  it  has  been  reported  that  data  justify  these  assumptions.  Accordingly,  R 
may  be  expressed  as 

n&<l>)  =  \^lT\<t>,p^)  (B4a) 

where  for  the  problem  of  plane  stress  (which  is  addressed  in  [B4]), 

and  r  is  magnitude  of  radius  vector  to  the  point  on  the  damage  surface  in  the  strain  space  de¬ 
fined  by  the  value  of  ^  and  two  surface  coordinates  Py  (7=  1,2).  The  subscript  1  refers  to  the  fi¬ 
ber  direction,  2  to  the  transverse  direction,  and  12  to  in-plane  shear. 

fi^=s^/\s\  ;  r  =  l,2  (B4c) 

In  [B4],  the  damage  surfaces  (in  the  three  dimensional  strain  space)  are  obtained  numerically 
from  analysis  of  data  from  notched  (±0)^^  specimens  for  various  values  of  0  and  different  loading 
paths  using  a  curve  fitting  procedure.  Programmed  and  controlled  loading  devices  measure  the 
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forces  as  well  as  the  displacements  at  the  loading  points,  so  that  the  total  energy  dissipation  in 
the  specimens  can  be  computed.  Finite  element  analysis  and  computed  elastic  strain  fields  in 
the  specimen  are  utilized  for  analysis  of  data. 

Authors  of  [B4]  plan  to  make  their  analysis  methodology  available  through  the  world  wide 
web,  but  currently  neither  the  method  nor  the  data  (damage  surfaces)  are  available  for  use.  For 
nonlinear  stress  analysis  of  laminated  plates  with  and  without  stress  raisers  (holes  or  cracks), 
the  following  simple  form  of  the  damage  surfaces  was  employed  in  [B5,  B6]. 

+ bi  i(<^)sgn  si}sf 

+{a22(?^)+b22(^)sgn^2}4  (B5a) 

+2ai2(^)^I^2 

=  1 

It  may  be  noted  that  (B5a)  is  a  Tsai-Wu  type  quadratic  failure  criterion  in  terms  of  strain 
(with  different  coefficients  in  the  four  quadrants  of  £1^2  plane  and  no  linear  terms).  The  coeffi¬ 
cients  (a,j  ±bn),  (a22±b22),ai2  and  333  are  functions  of  ^  and  they  may  vary  with  ^  in  different 
fashions,  but  subsequent  damage  surfaces  must  enclose  the  previous  ones.  Based  on  the  re¬ 
sults  of  studies  on  failure  theories  of  unidirectional  composites,  it  appears  that  the  representa¬ 
tion  of  the  damage  surfaces  by  (B5a)  may  not  be  very  accurate  for  small  values  of  <j>,  but  it  pro¬ 
vides  a  reasonable  estimate  of  stiffness  losses.  As  discussed  later,  the  representation  is  quite 
good  for  moderate  to  high  values  of  (j). 

The  residual  stiffnesses  of  Q-  after  some  energy  dissipation  obtained  using  the  constitutive 
law  are  given  below.  Note  that  Cy  in  (B2)  must  be  replaced  by  Qy  (reduced  stiffnesses  for 
plane  stress  problems). 

Qm  =Qii-(Aii(f^)+B„(«i)sgnfi) 

Q22  =Q22  -(A22(f^)  +  B22(f^)sgnf2)  (gg^.d) 

Qi2  =Qi2  ~Ai2(^) 

Q66  =Q66“A65(^) 

where 

[Ay(<#),By(^J)]  =  ^  [ay((^),by(<^)]d(#  (B5e) 

It  may  be  noted  that  (ay,by)  or  (Ay, By)  can  be  determined  if  stress-strain  responses  are 
known  from  displacement  controlled  tests  with  unidirectional  straining.  For  example,  Aee  (or  a^e) 
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will  be  known  from  the  nonlinear  shear  stress-strain  curve  for  the  unidirectional  composite. 
However,  unidirectional  straining  is  difficult  to  conduct  and  sudden  failures  may  often  occur  in 
testing  unidirectional  materials.  To  obtain  a  better  representation  of  the  in-situ  behavior  of  a 
lamina  in  a  laminate,  test  data  from  laminate  tests  are  preferable  for  determination  of  the  coeffi¬ 
cients  ajj  and  by  as  functions  of  (ji.  Tension  and  compression  test  data  of  (±0)ns  laminates 
(6=30, 45°, 60°)  were  utilized  in  [B5]  to  estimate  a22,  b22,  and  age  since  stiffness  losses  in  the  fiber 
direction  (reduction  in  Qn)  are  very  small  in  these  laminate  tests. 

Determination  of  au  and  bu  (or  axial  stiffness  losses)  is  a  more  difficult  problem,  because 
standard  tension  or  compression  tests  on  laminates  ((0/90)„3  or  (0/±45/90)  J  show  very  little 
nonlinear  response  up  to  the  maximum  load  at  which  sudden  failure  occurs.  The  Poisson  strain 
in  (0/90)„5  laminates  show  some  nonlinear  response  from  which  a^  was  estimated.  It  should  be 
noted  that  for  high  modulus  carbon/epoxy  composites  (such  as  AS4/3501-6),  the  transverse 
(extensional)  and  shear  stiffnesses  Q22  and  Qee  are  small  compared  to  the  axial  stiffness  Qn  and 
damage  initiation  strain  82  (or  yn)  is  usually  smaller  than  (or  of  the  same  order  as)  that  in  the 
fiber  direction  (si).  For  this  reason,  the  energy  to  be  dissipated  for  complete  loss  in  transverse 
(or  shear)  stiffness  is  much  smaller  (of  the  order  of  150-300  in  Ib/in")  than  that  required  to  cause 
complete  loss  in  fiber  direction  stiffness  (of  the  order  of  5000-7000  in  Ib/in®).  Displacement 
controlled  tests  on  (0/90),^  or  (0/±45/90)„3  laminates  are  very  difficult  to  conduct  and,  therefore, 
it  appears  that  the  best  approach  to  estimate  an,  bn  will  be  to  conduct  tests  on  laminates  with 
holes  or  notches,  and  then  correlate  the  results  with  those  from  appropriate  stress  (strain) 
analyses  using  various  chosen  values  of  an(^),bi,(^).  Matching  the  local  response  (load  vs. 
notch  opening  displacement  or  notch  tip  strain)  will  yield  the  desired  coefficients.  Calculations 
[B5]  show  that  it  is  possible  to  do  so  and  the  results  also  explain  the  so-called  “hole  size  effecf 
observed  in  tests.  It  may  be  noted  that  at  high  values  of  ^  (>300  in  Ib/in®),  all  the  coefficients  in 
(B5a)  except  an  and  b,,  and  all  the  stiffnesses  in  (B5b-d)  except  Qn  approach  zero,  but  the  loss 
in  fiber  direction  stiffness  continues  (initially  at  a  fast  rate  and  then  at  a  slower  rate)  up  to  very 
high  values  of  si  (>  5%),  which  indicates  that  a  lot  of  energy  is  dissipated  through  a  process  of 
progressive  fiber  breaks  (kinks)  followed  or  accompanied  by  gradual  fiber  pull  out  or  other  simi¬ 
lar  processes.  Therefore,  for  large  values  of  <l),  the  damage  surfaces  can  be  described  in  terms 
of  Si  alone.  In  many  problems  of  practical  interest,  it  may  not  be  necessary  to  model  the  mate¬ 
rial  response  up  to  such  high  values  of  dissipated  energy  density,  but  the  roles  played  by  inter¬ 
laminar  strains  and  stresses  may  be  quite  significant.  Responses  of  textile  composites  are  ob- 
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viously  influenced  by  such  stresses  and  strains  and  it  is,  therefore,  necessary  to  obtain  a  repre¬ 
sentation  of  including  the  effects  of  interlaminar  strains.  Such  representations  are  pre¬ 

sented  in  the  next  section.  It  should  be  noted  that  for  obvious  reasons,  no  influence  of  laminate 
stacking  order  is  predicted  by  the  model  discussed  above  [B5]  for  the  “hole  size  effect”  problem, 
and  the  representations  given  next  should  also  be  useful  for  estimating  the  influence  of  stacking 
in  such  problems. 
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DAMAGE  SURFACES  FOR  ATRANSVERSELY  ISOTROPIC  COMPOSITE 


The  problem  of  plane  stress  addressed  in  the  previous  section  is  an  approximation.  It  may 
be  noted  that  when  the  through  the  thickness  stress  0-3  is  zero,  it  is  highly  unlikely  that  the  cor¬ 
responding  strain  component  in  the  damaged  state  is  related  to  the  other  two  strains  and  £2 
via  the  elastic  relation,  i.e.,  the  sum  of  the  two  products  obtained  by  multiplying  these  strains 
with  the  appropriate  Poisson’s  ratios.  Such  a  relation  is  not  consistent  since  the  use  of  the  re¬ 
duced  stiffnesses  in  (B5b)  imply  that  when  0-2  is  zero,  £2  =  -Qu^i  which  may  not  neces¬ 
sarily  be  the  same  in  the  elastic  domain.  Experiences  in  the  use  of  the  theory  of  plasticity  indi¬ 
cates  that  inconsistencies  may  arise  if  a  plasticity  theory  is  formulated  by  choosing  some  of  the 
stress  components  as  zero.  Therefore,  it  appears  necessary  to  obtain  a  constitutive  law  retain¬ 
ing  all  the  strain  and  stress  components.  For  a  fiber  composite  which  is  transversely  isotropic  in 
the  elastic  state,  one  can  write  the  Helmholtz  free  energy  in  the  following  form. 

py/  +  2C^2£\£s  + 

+  GT(fd  +  X23)  +  GA(r?2+?'n)]  (B6a) 

where 


-  £’2  +  ^3 
=  ^2  “  ^3 


(B6b,c) 


k  is  the  plane  strain  bulk  modulus  and  Ga  and  Gt  are  the  axial  and  transverse  shear  moduli,  re¬ 
spectively.  The  stresses  are  given  by 

d{p\y) 

‘  ae, 


0-2+ 0-3  _d{pW) 
2  d£^ 

0-2  -  0-3  ^  ^PW) 
2  d£^ 


z-,2  = 

r,3  = 

^23 


^P¥) 

d(<P¥) 

^Yn 

d(P¥) 

5/23 


(B6d-i) 
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Based  on  the  assumption  of  a  quadratic  interaction  criterion  in  terms  of  stresses,  one  may 
choose  the  following  form  of  R  for  zero  (damage  initiation)  to  small  values  of  the  dissipated  en¬ 
ergy  density 


+a44(^)(ffd  +^23)‘‘‘®66(^)(x?2  +?'?3)] 


(B7a) 


where 


Si  -  Si  +ai2^s 


i?s=^s  +  «2i^i  (B7b-e) 

«12  =Cj2 

<^21  =  ^12  / 

Use  of  (B6a-i)  yields  a  transversely  isotropic  material  with  the  following  values  of  reduced  stiff¬ 
nesses  in  the  damaged  state  for  the  current  value  of  <j). 

Cj  1  =  C„  -  (A„  (4)  +  B, , (<^)sgn  £*i )  -  ( A22(<z>)  +  B22(?>)sgn  £*  )a^, 

Cl2  =Ci2  -(Aii(<^)  +  Bii(^)sgn4)«i2  -(A22W  +  B22Wsgn£s*)«2i 

k*  =  k-(A22(<^)  +  B22(<i>)sgnes)-(An(«>)  +  B,,((!))sgnf,)a:^ 

C44  =  Gj  —  A44(^) 

^66  =  Ga  “  A65(^) 


where 


[Aij(<^),Bij(^)]  =  (ay(<^),bij((^))d^>  (B7k) 

Representations,  which  are  more  complicated  than  (B7a),  may  have  to  be  chosen  to  match 
test  data  if  they  indicate  that  the  damaged  composite  behaves  as  an  orthotropic  or  a  more 
highly  anisotropic  (one  case  is  discussed  later)  material.  It  is  worthwhile,  however,  to  examine 
whether  the  simple  representation  can  yield  reasonable  correlation  for  small  values  of  (j).  Fur¬ 
ther  simplification  can  be  achieved  by  choosing 

b22(<^)  =  a22(<^) 

B22(^)  =  A22(<^) 

which  imply  that  hydrostatic  compression  in  the  plane  of  transverse  isotropy  does  not  cause  any 
damage,  or  in  other  words,  damage  may  occur  only  at  very  high  values  of  hydrostatic  compres¬ 
sive  strain,  which  is  not  of  interest  in  problems  of  the  type  under  consideration  [B7,  B8].  Fur¬ 
ther,  for  small  to  moderate  values  of  <[),  it  is  possible  to  use  the  following  representations  for  the 


(B7l,m) 


losses  in  various  stiffnesses  in  (B7f-k)  and  (B7m)  to  match  test  data  for  softening  observed  in 
(±0)„  and  other  layups. 


2 A22(^)  =  A22(<^)  +  B22((^)  =  k  1  -  (l  +  ^  /  ^>2  )  "' 

L  J  (B8a-d) 

A44((#)  =  GT[l-(l  +  f^/(^4p] 

A66(^>)  =  Ga[i-(1  +  ^»/^6P 

It  may  be  noted  that  these  representations  are  adequate  when  the  stress-strain  relations  are 
of  the  power  law  type  after  the  onset  of  damage.  For  example,  if  the  axial  shear  stress-strain 
relation  is  of  the  following  form  {tx2,y\2  being  the  stress  and  strain  at  damage  onset) 


^12  “  ^?2(?'l2  /7?2) 


(B8e) 


then  n6  and  (t)6  are  given  by 
06  =2/(^6 +1)-1 
=GAr?206/2 

From  test  data  for  AS4/3501-6  laminates  reported  in  literature  and  fits  performed  in  [B5] 
(where  ajj(<i5),bjj((z>),Aij(^),Bij(<J),  etc.,  were  chosen  by  linear  interpolation  in  (j)  from  tabulated 

values),  it  appears  that  the  following  parameters  should  yield  reasonable  fits  to  the  laminate  re¬ 
sponses  when  <p  is  not  iarge. 

nf- 2.576,  <z)^'=4812inlb/in^ 

nr=  1.372,  ^f=2810inlb/in^ 

nj  =  2.538,  =22.68  in  lb /in^  (B9a-e) 

n4  =  0.287,  =11.40inlb/in^ 

n6=  0.478,  =  14.86  in  lb /in^ 

Representations  of  the  types  described  in  the  sets  of  equations  (B7,  A8)  are  currentiy  being 
used  for  problems  with  small  to  moderate  stiffness  losses  in  some  textile  composites. 

Power  law  type  representations  with  constant  values  of  the  indices  nj  may  not  be  always 
adequate,  especially  for  larger  values  of  ^ .  Better  representations  can,  however,  be  obtained 
by  choosing  the  indices  ni  as  functions  of  ^ .  For  example,  one  may  choose 

A„  W)  ±  B„(«  =  C„  []  -  (l  -  ^  /  ,>f  r" ’1  (B1 0) 


A„(<^)±B„(<^)  =  C„  \-(\-(l>l4>\) 
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instead  of  (B9a)  and  similar  relations  for  the  quantities  in  (B9b-d).  Differentiation  of  equations  of 
the  type  (B10)  yields  the  coefficients  aij(^),bij(^).  For  example,  from  (B10)  one  can  obtain 

where 

„±'  dnf  (B12) 

Although  representations  of  the  type  discussed  here  should  be  adequate,  when  n;  are  weak 
functions  of  ^ ,  it  is  necessary  to  check  that  the  damage  surfaces  for  any  value  of  enclose 
the  ones  for  values  of  . 

It  was  pointed  out  in  the  last  section  that  for  large  values  of  ^ ,  the  damage  surfaces  can  be 
described  by  alone,  since  matrix  dominated  stiffnesses  reduce  to  negligible  values  for  mod¬ 
erate  values  of  (j) ,  but  fiber  direction  stiffnesses  reduce  at  a  slower  rate  as  straining  continues. 
Because  of  this  behavior,  one  should  expect  that  Uxi  chosen  in  (B7a,b)  as  a  constant  should 
change  with  ^  reducing  to  zero  for  moderately  large  values  of  ^ .  One  should  note,  however, 
that  values  of  o:,2  =C,2  /C„  for  high  modulus  fiber  composites  are  usually  small  and  hence  a 
modified  and  integrated  version  of  (B7a)  given  below  may  be  better  suited  for  modeling  the  be¬ 
havior  almost  to  the  point  of  complete  stiffness  loss.  The  integrated  form  is  useful  when  ni  are 
functions  of  ^  (as  in  (B10-A12)). 

R(£,  (f)A(l)  =  ^[(A,  ,(^)  +  B,  1  sgn  f ,  )e\ 

+  (A22(<z5)  +  B22Sgnf*)£*^  {B13a) 

+A44(^)(fd  +  X23)  +  A66(^)(^?2  +  >'b)] 

It  may  be  noted  that  if  is  replaced  by  s\  and  (B13a)  is  differentiated  with  respect  to  ^  one 
obtains  (B7a).  Since  a, 2  is  small,  results  obtained  using  (B13a)  should  not  differ  significantly 
from  those  using  (B7a)  for  small  stiffness  losses.  It  should  also  be  noted  that  the  choices  (B7I, 
m)  will  not  be  suitable  for  large  values  of  ^ ,  since  stiffness  loss  would  be  expected  under  in¬ 
plane  hydrostatic  compression  at  large  strain  levels.  Therefore,  one  additional  coefficient  must 
be  chosen  by  matching  data  from  tests  under  appropriate  strain  states/paths.  Use  of  (B13a) 
and  (B6a-i)  yields  reduced  stiffnesses  in  the  forms  (B7f-j)  with  the  following  changes  in  (B7f-h). 
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(B13b-d) 


Cl  1  =  C 1,  -  ( A„  (<^)  +  B,  1  (^) sgti  )  -  ( A22  +  B22  sgn  )a2 1 
C*i2  =  C12  -(A22  +  B22  sgnf*)a2i 

k*  =  k-(A22+B22Sgnf*) 

If  one  wishes  to  simulate  states  approaching  complete  stiffness  loss  (for  large  values  of  ^),  one 
may  choose  the  following  representations  instead  of  (B8a,b). 

An((#)±Bi,(^)  =  (C„-Ci2/k)  l-(l  +  ^5/$5^') 

L  J  (B13e,f) 

A22(«>)±B22(<^)  =  + 

A44((^),A66(<!>)  can  be  chosen  to  be  the  same  as  (B8c,d).  The  indices  nf  may  be  assumed  to  be 
weak  functions  of  (j)  as  discussed  earlier. 

Finally,  one  should  note  that  the  representation  (B13a)  involving  squares  of  the  strain  vari¬ 
ables  and  those  of  the  types  (B10-A12)  are  still  very  simple  and  may  not  be  adequate  in  some 
cases.  One  phenomenon,  which  may  need  some  consideration,  is  the  effect  of  in-plane  hydro¬ 
static  strain  (or  stress)  on  damage  in  a  state  dominated  by  shear.  Experimental  data  indicate 
that  shear  strains  (or  stresses)  to  cause  initial  (and  progressive)  damage  are  usually  increased 
when  hydrostatic  compression  [B7,  B8]  is  present.  Such  behavior  can  possibly  be  modeled  by 
choosing  the  following  modified  form  of  (B13a). 

\t  Rfe  =  ^[( Ai  1  i<l>)  +  Bi  1  Sgn  )sl 

+  (A22(<^)  +  B22Sgn£sK^ 

+  A44(^)(ed 23)(1 

+A66(^)(ff2  +7b)(1  +  ^2(^)^®S^2)] 

where 


t^nO,=(sl-,rhrfs:  (B14b) 

tan  02  =  (rn+rh 

and  5ii^),S2(.^)  functions  of  ^  to  be  chosen  to  match  test  data.  It  is  expected  that  Si  and 
S2  should  approach  zero  for  large  ^  and  accordingly  they  may  be  chosen  as 

#2(«>)  =  'S2o(l  +  <S/|*6r"‘“ 
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^10  and  S20  may  be  selected  based  on  available  data  for  initial  damage  following  procedures 
similar  to  those  used  in  [B8].  07  and  ns  may  be  chosen  if  subsequent  damage  surfaces  under 
combined  shear  and  hydrostatic  in-plane  compression  and  tension  can  be  estimated  from  test 
data.  As  initial  estimates  07  and  ng  may  be  chosen  equal  to  tu  and  ne,  respectively.  It  may  be 
noted  that  evaluation  of  stresses  using  the  representations  (B14a-e)  and  the  relations  (B6d-i) 
indicates  that  the  stress  total  strain  relations  in  the  damaged  state  are  complicated  and  the 
stiffness  matrix  (reduced)  is  no  longer  symmetric.  This  is  a  resuit  of  the  complex  but  more  re¬ 
alistic  damage  surfaces  chosen  (B14a)  and  has  been  discussed  in  [B4].  This  complexity  also 
necessitates  the  examination  of  the  damage  surfaces,  so  as  to  ensure  that  the  present  one  en¬ 
closes  the  previous  surfaces. 
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DAMAGE  SURFACES  FOR  AN  ISOTROPIC  MATERIAL 


In  addition  to  the  unidirectional  fiber  composites  in  wavy  forms,  textile  composites  contain 
matrix  pockets,  which  are  also  brittle  and  will  undergo  stiffness  loss  with  increasing  strain.  The 
Helmholtz  free  energy  for  such  materials  can  be  expressed  in  the  following  form  based  on  the 
principles  of  damage  mechanics  and  assumptions  similar  to  those  used  for  the  damage  sur¬ 
faces  (in  the  previous  section)  for  fiber  composites. 
pil/=\  [k4  +Gyefr]-^  R(£,^)d<^ 

/eff  =  2(fn  +  +  £33)  +  f  +  ?'23  + 


(B15a,b) 


where  K  is  the  bulk  modulus,  G  is  the  shear  modulus, 
f V  =  1  +  ^22  +  ^33»  volumetric  strain 

fjj  =  fjj  -  fy  /  3;(i  =  j),  deviatoric  strains  (B1 5c-e) 

7ij  =  shear  strains;  (i  ^  j) 

R(f,^)  =  1  (for  damage  growth  to  occur)  (B1 5f) 


(B15f)  define  the  damage  surfaces  for  various  values  of  ^  and  no  damage  growth  will  occur  if 


R(^^)<1. 

The  stresses  are  given  by 

o’m  =  (o’!  1  +  o’22  +  o’33)  / 3  =  niean  stress 

dSy 

cr|:  =  an  -  cr„  =  — (i  =  j),  deviatoric  stresses 

ij  u  m 

J-  =  .  (j  ^  shear  stresses 

dr-ii 

The  simplest  form  of  R(f,^),  which  may  be  chosen  is  of  the  following  forms 


(B15g-i) 


=  ^[(av(«^)  +  by(<^)sgn  fy  )4  +  as(<^)reff 
Ri£,^)d^  =  ^[(Ay(^))+ By(^>)sgn  fy)4  +  As(<^)/eff] 


(B16a,b) 


Substitution  of  (B16b)  in  (B15a)  and  use  of  (B15g-i)  yield  the  following  expressions  for  stresses 

c7ij=|K-(Ay((i)  +  By((!5)sgnfy)-|(G-A,((;)))|fy^ij 

+  2(G-A3(^^))£ij;i=j  (B16c,d) 

^ij  “  {G  “  ij’  *  ^  j 
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It  is  clear  that  the  damaged  material  behaves  as  an  isotropic  material  with  a  reduced  bulk 
modulus  of  {K-(Av(<i))±Bv(<^)}  and  a  reduced  shear  modulus  of  {G-As(^)}.  More  compli¬ 
cated  representations  can  be  chosen,  if  desired,  for  the  purpose  of  matching  test  data  for  vari¬ 
ous  strain  paths  and  strain  ranges,  but  it  appears  worthwhiie  to  examine  the  usefuiness  and 
range  of  validity  of  the  simple  representations.  Additionai  simplifying  assumption  may  be  made 
by  neglecting  the  effect  of  hydrostatic  compression  (vaiid  over  a  limited  range  of  strains)  and 
choosing 

=  (B16e.f) 

B,(<#)  =  A,(^) 

Based  on  power  law  type  representation  of  stress-strain  relations  (discussed  for  composites  in 
the  previous  section)  beyond  the  elastic  region,  one  may  also  choose 


(B17a,b) 


2A^(«))  =  =  Kp  +  ({^  / 

A3(«))  =  g[i-{1  +  ^>/<^3)""‘] 


with  Hv  and  nj  chosen  as  constants  (for  low  to  moderate  values  of  ^)  or  as  weak  functions  of  ^ 
(for  wider  range  of  strains  and  ^ ).  In  the  latter  case,  however,  effect  of  hydrostatic  compression 
may  have  to  be  considered  (assumptions  (B16e,f)  may  not  be  valid)  and  it  would  be  necessary 
to  introduce  additional  parameters,  i.e., 

Av(<^)±B^(«J)  =  K  I  (B18) 


Finally,  the  increase  in  shear  strains  (or  stresses)  to  cause  damage  with  increase  in  hydro¬ 
static  compressive  strain  (or  stresses)  can  be  modeled  by  the  following  modification 

t  =  rflAvCf^) + Bv(<^)sgn 

21^  >  (B19) 

+  As(<^)reff(l  +  ^vCOS^)] 

where 

tan^  =  -^  (B20) 

and  (a  constant  or  a  function  of  ^)  has  to  be  chosen  to  match  experimental  results  similar  to 
those  used  to  select  the  parameters  in  Mohr-Coulomb  type  failure  theories  (see,  for  example, 
[B8]).  The  resulting  stress-total  strain  relations  are  complicated  and  the  stiffness  matrix  may  not 
be  symmetric  as  discussed  in  the  previous  section. 
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APPENDIX  C  -  THERMO-MECHANICAL  PROPERTIES 
OF  A  REPRESENTATIVE  AREA  ELEMENT 


SUMMARY 

This  Appendix  examines  fundamental  micromechanical  relations,  and  how  they  may  be 
generalized  to  cover  a  wider  class  of  material  forms.  In  particular,  the  typical  homogenization 
process  applied  to  multiphase  materials  has  been  modified  into  a  form  more  amenable  to  textile 
composites.  This  approach  supports  a  goal  of  obtaining  useful  results  using  analytic,  closed- 
form  methods.  The  methodology  can  also  be  supported  by  finite  element,  and  other  numerical 
methods.  But  the  numerical  methods  must  be  supported  by  strong,  fundamental  mechanics. 

INTRODUCTION 

Laminated  composite  plates  containing  layers  of  unidirectionally  reinforced  continuous  fiber 
composites  of  various  orientations  are  now  being  widely  used  for  various  structural  applications. 
The  state-of-the-art  of  evaluating  their  structural  performance  is  quite  advanced  and  their  me¬ 
chanical  behavior  is  well  defined.  On  the  other  hand,  the  mechanics  of  woven  fabric  reinforced 
composite  plates,  manufactured  by  combining  the  age  old  technology  of  textiles  and  processing 
methods  recently  developed  for  laminated  fiber  composite  plates,  are  still  not  clearly  under¬ 
stood.  Several  attempts  have  been  made  to  characterize  the  structural  response  of  such  plates 
based  on  various  geometric  parameters  and  properties  of  the  constituent  fiber  and  matrix  mate¬ 
rials  (for  a  survey,  see  [C1]).  The  main  idea  behind  these  studies  is  to  use  the  concept  of 
equivalent  homogeneity,  widely  used  in  mechanics  of  fiber  reinforced  composites  [C2].  These 
methods  attempt  to  find  the  relationship  between  either  average  stresses  and  average  strains 
or  average  stress  resultants  and  average  mid-plane  strains  and  plate  curvatures.  We  first  con¬ 
sider  the  definitions  of  these  average  quantities  and  give  some  simple  theorems  relating  them 
to  displacements  and  tractions.  The  displacements  and/or  tractions  are  applied  on  the  bound¬ 
ary  of  a  representative  area  element  (Figure  C1),  which  consists  of  several  nonhomogeneous 
layers.  Each  layer  is  either  a  homogeneous  material  (such  matrix  layers  are  often  used  in 
practice)  or  a  fabric  reinforced  composite.  Each  fabric  reinforced  layer  can  be  further  divided 
into  a  number  of  sublayers,  if  needed.  It  is  assumed  that  each  layer  (or  sublayer)  consists  of  a 
finite  number  of  materials  (unidirectionally  or  bidirectionally  reinforced  composite  or  matrix  ma¬ 
terial),  whose  thermomechanical  properties  are  known.  For  calculating  the  thermomechanical 
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Figure  C1 .  Representative  Area  Element  (RAE)  Containing  Several  Layers 
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DEFINITIONS  AND  THEOREMS 


Average  Displacements  and  Forces 

We  define  the  average  mid-plane  strains,  curvatures  and  transverse  shear  strains  in  the 
following  manner: 


h/2 


£::  = - 

2hA 


Ja  j 


-h/2 

h/2 


=  U3 


I  h/2 


;  ij  =  1.2 


(Cl) 


ri3=2%=— Ja  J  (ui,3  +  U3,i)dzdA 

hA  _h/2 


where  h  is  thickness  of  the  plate  and  A  is  the  area  of  representative  element  (RAE).  An  exam¬ 
ple  is  shown  in  Figure  C1.  ui,  U2  are  the  displacements  in  x  and  y  directions  which,  in  general, 
depend  on  the  space  variables  x,  y  and  z.  Obviously,  the  size  of  the  RAE  should  be  large  com¬ 
pared  to  the  “mini”  structural  dimensions  to  make  the  averaging  process  meaningful  [C1].  If 
these  quantities  are  to  be  used  for  determining  the  plate  response  in  practical  applications,  the 
size  of  the  structural  plate  element  should  be  sufficiently  large  as  compared  to  that  of  the  RAE. 
Further,  the  averages  are  meaningful  when  they  do  not  fluctuate  with  increasing  size  of  the 
RAE. 

By  the  application  of  divergence  theorem  and  continuity  conditions  for  displacements  at  the 
interfaces,  it  can  be  shown  that 


-  1  '  /  \ 

J  (uinj  +  UjiijjdzdC 


-h/2 

h/2 


U  3 


Jc  j  (uinj  +  Ujni)zdzdC 


h^A'''  -m 


:  iJ  =  1.2 


(02) 


h/2 


j  U3nidzdC  +  j  {uj(h/2)-Ui(-h/2)}dA 
-h/2  A 


where  tij,  n2  are  the  direction  cosines  of  the  normal  at  a  point  on  the  boundary  C  of  the  RAE. 
Average  in-plane  stress  resultants,  moments  and  shear  forces  are  defined  as; 
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__  I  h/2 

Nij=-7l  ^  o-jjdzdA 

A  A  -h/2 

_  I  h/2 

My  =  —  J  j  cTjj  zdzdA 

A  A  -h/2 
__  1  h/2 

Qij=-r^  J  o-jsdzdA 

A  A  -h/2 


:  ij  =  1 .2 


By  the  use  of  the  equilibrium  equations,  the  divergence  theorem  and  the  continuity  condi¬ 
tions  for  tractions  across  interfaces,  the  average  stress  resultants  and  moments  can  be  related 
to  the  boundary  tractions  on  the  boundary  C  and  top  and  bottom  surfaces  (z  =  ±h/2)  of  the  plate. 
These  relations  in  symmetricized  forms  are  given  below: 

f  T(xjTi+XiTj)dzdC 
^A[_c  -h/2 

+  J  {x,[T,(+h/2)  +  T,(-h/2)]+x,[Tj(+h/2)  +  Ti(-h/2)]}<IA  ] 

A 

1  r/  \  1 

*^0=7^^  f  {  xjTj+XjTj  z-XjXjTjjdzdC 

2A|-C-h/2‘  .ij  =  1,2  (C4) 

/  {^[Ti(+h  /  2)  -  T;(-li  /  2)] + i^[T,(+h  /  2)  -  Tj(-h  /  2)] 

-x,Xj[T3(+h/2)+T3(-h/2)]]}dA 

_  1  r  h/2  ,  - 

Qi  =—  J  j  XjTj  dzdC+  J  Xi[T3(+h/2)+T3{-h/2)]dA 
ALc  -h/2  A 

It  may  be  noted  that  when  the  surfaces  z  =  ±h/2  are  stress  free,  the  second  integral  over  the 
area  A  on  the  right  hand  side  of  each  of  the  above  set  of  equations  as  zero  and  the  average 
values  NjjjMjj  and  Qj  are  expressed  in  terms  of  integrals  evaluated  over  the  boundary  C.  The 

same  result  is  approached  for  sufficiently  large  RAE  when  the  surface  tractions  are  oscillatory 
and  periodic  functions  of  in-plane  coordinates  Xj,  X2.  In  addition,  if 


h/2 

J  J  XiXjT3dz  =  0 

C  -h/2 


(C4a) 


then  Njj,M  y  are  obtained  from  T|  and  T2  on  boundary  C. 


Constitutive  Relations 


For  fabric  reinforced  plates  containing  repeating  arrangement  of  weave  patterns  of  the  type 
shown  in  Figure  C1,  it  is  expected  that  the  average  stress  resultants  and  moments  will  be  re- 
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lated  to  the  mid-plane  strains  and  curvatures  by  constitutive  relations  similar  to  the  ones  used  in 
laminated  plate  theory.  For  linear  elastic  problems,  these  relations  are  of  the  form  given  below. 


N,' 

'Nn' 

All  Ai2  Ai6 

N2 

N22 

A 12  A22  A26 

N,2 

A 16  A26  A56 

M, 

M„ 

M2 

M22 

Sym. 

M6_ 

_Mi2. 

where 


^1 

^11 

^2 

^22 

^1 

2£|2 

^11 

K2 

^22 

.^6. 

2k  12  _ 

(C5a) 


and  the  last  vector  on  the  right  hand  side  of  (C5a)  quantifies  the  stress  resultants  and  moments 
due  to  change  in  temperature  in  absence  of  strains  and  curvatures.  The  relations  (C5a)  are 
obtained  based  on  the  following  in-plane  displacement  fields 

Ui=(fij  +  z^rij)xj  (C5b) 

and  the  interlaminar  stress  fields 

O-zz  =  ^X2  =  ^■yz  =  0 

in  each  layer  of  the  laminate,  which  yields 


(A;,B*j,Dy)=  J  C(j(l,z,z2)dz 

— h/2 

(ri,Ai)=  j /ij(l,z,z)ATdz  ;ij  =  1,2,6  (C5d) 

-h/2 

Cy  =  Cy  -  Ci3Cj3  /  C33 

where  the  constitutive  law  for  each  lamina  is  of  the  type 

o-j  =Cijfj+riAT  :ij  =  1,2,  ...6  (C5e) 
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The  strain  energy  in  area  A  of  the  plate  can  be  written  in  the  following  form  in  contracted  nota- 


+2(ri  fj  +  Ai/fjj-Frj 
Ft=  T((r3ATf /C33)dz 


Midplane  strains  and  curvature  due  to  a  temperature  change  aione  are  expressed  in  terms  of 
thermal  expansion  and  curvatures  given  by: 


«6  A  B  Tg 

I3\  [B*  D*J  a’i 

A*  and  D*  matrices  in  equation  (C5)  and  (C6)  are  the  effective  in-plane  stiffnesses  and  bending 
rigidities  and  the  B*  matrix  contains  the  bending-extension  coupling  terms. 

The  shear  forces  are  related  to  the  average  shear  strains  by  relations  of  the  following  type 

Qi  _  ^11  Ym 

Q2  Ki2  K22  723. 

In  subsequent  stages  of  the  work,  we  will  study  the  effects  of  damages  (crack-like  defect) 
and  their  growth  due  to  increasing  load.  In  such  cases  the  elements  of  the  matrices  A*,  B*  and 
D*  may  be  assumed  to  depend  on  the  state  of  strain  in  the  RAE. 


EVALUATION  OF  EFFECTIVE  ELASTIC  PROPERTIES 


To  evaluate  the  effective  elastic  properties,  the  standard  technique  is  to  first  apply  a  homo¬ 
geneous  boundary  condition  in  terms  of  displacements  or  tractions  [C2]  on  the  heterogeneous 
element  (RAE).  Sometimes  it  is  advantageous  to  apply  mixed  conditions,  i.e.,  some  in  terms  of 
displacements  and  some  in  terms  of  tractions.  If  exact  solutions  can  be  found,  then  the  strain 
or  complementary  energies  in  the  element  are  equated  to  that  in  a  homogeneous  material  to 
calculate  the  effective  elastic  properties.  However,  the  geometric  arrangements  of  the  con¬ 
stituents  are  often  so  complicated  that  exact  solutions  are  difficuit  to  obtain  and  approximate 
solutions  (simple  or  finite  element  calculations)  are  usually  employed.  In  what  follows  we  dis- 
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cuss  some  boundary  conditions  and  give  some  approximate  solutions  which  can  yield  bounds 
on  effective  properties  as  well  as  reasonable  estimates  of  stress  and  strain  fields  in  the  RAE. 


Boundary  Conditions 

We  consider  the  following  alternative  types  of  boundary  conditions  on  the  surface  S  (con¬ 
sisting  of  C  and  z  =  ±h/2)  of  the  RAE.  These  choices  are  motivated  by  the  displacement  fields 
in  a  laminated  plate  due  to  prescribed  midplane  strains  and/or  curvatures. 


Set1 

;i  =  l,2  onC 

on  C  (C8a) 

on  z  =  ±h  /  2 

Set  2 


T3=0 

Ti=T2=T3=0 


Ui=(fg  +  2rH)xj  ;i  =  l,2 

T3  =  0  on  S 


(C8b) 


Set  3 

Ui={fP  +  z?f|)xj  ;i  =  l,2 

j  (C8c) 

U3  =  -2^uXiXj  +  F(z)  onS 

where  F(z)  (constants  at  2±h/2  and  a  function  z  on  C)  is  yet  to  be  chosen.  It  may  be  noted  that 
use  of  (C2)  and  the  third  of  equation  (C4)  yields 

Ki^=Kl  .  .  _^  2  (CQa.b) 

Qij=o 


for  the  boundary  conditions  (C8a)  or  (C8b),  whereas  for  the  case  (C8c) 


/i3=0 


(C9c) 


Therefore,  they  can  be  used  to  obtain  the  effective  properties  which  relate  the  average  stress 
resultants  and  moments  to  average  midplane  strains  and  curvatures. 

For  AT=0,  the  strain  energy  in  the  RAE  can  be  written  in  the  following  form 
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When  AT?iO,  the  following  must  be  added  to  (C10) 

+T3^''u3)dS  (011) 

The  integrals  in  (C10)  and  (C11)  are  evaluated  on  surfaces  where  the  tractions  are  not  pre¬ 
scribed  as  zero  (as  in  (C8a)  or  (C8b)).  Using  (C8a-c),  (C4),  (C10)  and  (C1 1)  one  obtains 

U^=A[NijfP  +  M..^P+C,]  forAT=0  (C12a) 

=  A[Nf  +  MfVg  +  Cf"^]  for  AT^^O  (C1 2b) 

where 

Cj=o  for  cases  (C8a,b)  (C13a) 

h/2 

C,  =  J  J  F'(z)cr33dzdA/A  for  case  (C8c)  (C13b) 

A  -h/2 

Ny,  M  y  are  related  to  the  tractions  by  (C4).  The  quantities  with  superscript  AT  are  due  to  tem¬ 
perature  change  alone.  If  the  woven  or  fabric  reinforced  plate  has  a  constitutive  law  similar  to 
(C5),  then  (C12)  and  (C13)  indicates  that  the  strain  energy  is  of  a  form  similar  to  (C5f).  It  may 
be  noted  that  the  term  Ci  or  cf^  for  case  (C8c)  appear  because  of  the  integral  jT3F(z)ds  and 
Ts^iO  in  this  case.  The  last  of  (C13)  is  obtained  by  the  use  of  the  divergence  theorem  and  the 
equilibrium  equation  in  terms  of  stresses  (in  z-direction)  in  the  following  manner 

,  h/2 

jT3U3dS  =  J  jT3F(z)dzdC+|  (T3(h/2)F(h/2)  +  T3(-h/2)F(-h/2))dA  (C14a) 

C-h/2  A 


and 

J  JT3F(z)  dC  dz  =  J  J  r3jt;3F(z)  dC  dz 
=  .f^^3i,iF(z)dzdA 

=  -  J  J  0-33  3F(z)  dz  dA 

h/2 

=  -  J  (T3(h  /  2)F(h  /  2)  +  T3(-h  /  2)F(-h  /  2))dA  +  J  J  F'(z)c733  dz  dA 

A  A  -h/2 

The  first  term  on  the  last  expression  of  (C14b)  cancels  with  the  second  term  in  (C14a)  and  we 
obtain  the  expression  Ci  as  it  appears  in  (C13b). 


A  Simple  Upper  Bound 


For  a  displacement  field,  we  assume  that  equation  (C8c)  also  holds  inside  the  RAE  and  F(z) 
is  yet  an  unknown  function  of  z  and  its  derivative  is 

f(z)  =  dF(z)/dz  (C15) 

It  follows,  therefore,  that  the  state  of  strain  in  the  RAE  is  given  by 
£jj  =  £P+z/trg  ;iJ  =  U2 

f33  =  f(z)  (C16) 

and 

;i  =  l,2 

and  the  totai  potential  energy  for  the  assumed  field  can  be  written  in  the  foliowing  form  using 
the  contracted  notation 

T  [Cy(4^j^+2z4/f?+/rP/c?) 

^  A-h/2 

+  C33f^(z)  +  2Cj3(£:f +  z/(rf)f(z)  (d7) 

+  |27'j(£'f  +  z/cP)  +  273f(z)|AT]dz  dA 

where  the  Qj  and  y,,  are  the  thermoelastic  constants  which  are  functions  of  position  and  their 
values  depend  on  the  constituent  r  with  the  following  constitutive  iaw  (similar  to  C5f) 

<Ti  =  Cij£j +7iAT  i>j  =  1.2,3,6  (C18) 

It  may  be  noted  that  for  a  laminated  plate,  equation  (C17)  reduces  to  the  energy  expression  as 
per  the  laminated  plate  theory  solution  provided  f(z)  is  in  each  lamina  is  chosen  so  as  to  make 
033  =  0  in  that  lamina  for  the  strain  state  (C16).  For  the  plate  with  a  woven  construction,  we  note 
that 

U®>U  (C18) 

where  If  is  given  by  (C17)  and  U  is  the  actual  strain  energy  (or  potential  energy)  which  must  be 
a  minimum.  We  can  now  write  (C17)  as 

^  -h/2^ 

+  C33f^(z)  +  2Cf3(£j+z/rf)f(z)  (C19) 

+  2?^?  (fP  +  z/rf)  +  2;'3f(z)]  dz 

where 
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(C20) 


cfj  =  J  CydA/A 

A 

rf  =  j  ^'jATdA/A 

A 

We  may  now  seek  a  minimum  value  of  If  for  determining  f(z)  which  yields 


f(z)  =  - 


+  ZKf)  +  r3 


which  also  implies  that 

J(T33dA  =  0  for  all  z 

and  therefore  from  (C13) 

for  case  C8c. 

Substitution  of  (C21)  in  (C19)  yields  the  following  expression  for  U® 

+  2(r/4H-A'i^P)-Fj  ] 


(C21) 


(C22) 


(C23) 


where 

h/2 

(A|j,Bij,Dy)=  I  Cy(l,z,z^)dz 
“h/2 

(ri',Ai)=  J  y\{\,z)Az 

-h/2 

cjj'  =  Cfj  -  Cf3C?3  /  CI3  :  ij  =  1 ,2,6  (C26) 

Y\=y■.-^^y%lch 

and 

Ft=  jJ(r3f 

These  results  yield  upper  bound  type  estimates  of  the  desired  effective  thermoelastic  constants 
for  the  plate  and  they  can  be  evaluated  for  any  known  geometric  arrangement  obtained  by 
stacking  fabric  reinforced  layers  with  or  without  stitching  or  other  through  the  thickness  rein¬ 
forcements.  It  may  be  noted  that  for  the  case  when  several  fabric  reinforced  layers  are  stacked 
up,  it  is  only  necessary  to  evaluate  the  effective  constants  for  one  layer  and  then  using  the  fol¬ 
lowing  relations  to  compute  the  properties  for  the  plate. 


C-10 


[A]  =  Z  [a"'] 

[B]  =  i:  [a"']z^+[b"^] 

[D]  =  z  [a"'](z^)''  +2[b"^]z^  +[d"^] 

Where  the  summation  is  taken  over  all  layers  (denoted  with  superscript  k)  and  is  the  dis¬ 
tance  of  the  midplane  of  the  layer  k  from  the  midplane  of  the  plate.  It  may  be  noted  that  the 
results  given  above  also  yield  upper  bound  estimates  for  the  whole  plate  since  the  displacement 
compatibility  is  satisfied  everywhere.  These  results  can  be  viewed  as  a  generalization  of  results 
obtained  for  effective  moduli  of  fabric  layers  based  on  the  assumption  of  iso-strain  conditions 
[C1]. 

More  accurate  upper  bound  estimates  can  be  obtained  by  choosing  other  types  of  dis¬ 
placement  fields  in  the  RAE  subjected  to  boundary  conditions  (C8a),  (C8b),  or  (C8c).  Condition 
(C8a)  imposes  less  restrictive  conditions  on  displacements  and,  therefore,  may  be  better  suited 
to  allow  simple  choices  of  other  fields.  However,  since  the  geometry  is  usually  quite  compli¬ 
cated,  it  may  be  necessary  to  obtain  numerical  solutions  for  nonzero  value  of  one  prescribed 
quantity  such  as  s\  and  zero  values  for  others  and  evaluate  the  total  strain  energy  in  all  the 
constituent  regions  (or  by  evaluating  Nij  and  Mij  by  using  (C4)). 
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APPENDIX  D  USE  OF  P-ELEMENTS  FOR  STRESS  ANALYSIS  OF  A  TEXTILE 


INTRODUCTION  AND  SUMMARY 

Additional  innovative  analysis  approaches  were  also  investigated.  One  approach  was  to 
use  P-element  technology  to  model  textiles.  P-elements  use  increasing  polynomial  order  to 
obtain  convergence.  A  major  advantage  of  the  approach  is  that  single  elements  can  be  used 
to  represent  large  portions  of  the  textile  architecture,  potentially  simplifying  the  modeling.  A 
commercial  P-element  code  called  Mechanica  [D1]  was  used  to  test  the  approach.  This  pro¬ 
vides  a  high  resolution  picture  of  the  linear  stresses  within  the  textile  unit-cell.  However,  the 
modeling  tools  available  in  the  commercial  package,  while  more  expedient  than  tradition  FE 
meshing  methods,  were  too  time  consuming,  and  are  difficult  to  automate.  This  drawback, 
combined  with  the  limitation  that  the  P-element  method  cannot  handle  general  nonlinear  mate¬ 
rial  properties  has  lead  us  to  focus  on  tradition  FE  methods. 

DISCUSSION 

The  P-element  model  is  based  on  the  same  geometric  assumptions  used  in  the  automated 
meshing  routines  described  above.  The  geometry  was  transferred  to  Mechanica  as  a  series  of 
curves,  and  the  mesh  created  using  the  graphical  tools  provided  in  the  software.  Figures  D1 
and  D2  show  portions  of  the  resulting  model.  The  individual  P-elements  can  form  into  the 
shape  of  the  undulating  section,  allowing  for  a  comparatively  coarse  mesh.  General  repeating 
boundary  conditions  are  not  possible  in  this  code  (due  to  the  lack  of  a  mult-point  constraint  ca¬ 
pability).  By  using  symmetry  conditions,  the  axial  load  condition  can  be  applied. 

Stress  results  for  a  segment  of  warp  yarn  are  shown  in  Figures  D4-D6,  using  the  coordi¬ 
nate  system  shown  in  Figure  D3.  The  results  are  for  a  fixed  displacement  in  the  warp  direction, 
resulting  in  an  average  on  stress  of  22.9  GPa.  These  figures  show  how  the  undulations  give 
rise  to  033,  and  T13  stresses  which  can  cause  intra-bundle  cracking.  In  addition,  there  is  a  con¬ 
centration  of  On  stress  that  will  reduce  the  fiber  failure  loads  as  compared  to  a  conventional 
laminate.  The  maximum  stresses  in  a  warp  yarn,  normalized  by  the  unit  cell  average  stress,  are 
summarized  in  Table  D1 . 
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Figure  D1 .  P-Element  Model  of  Yarns  for  Plain-Weave  Unit  Cell 
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D-3 


Stress  XZ  {Top) 

Avg.  Max  +2.0538E+00 
Avg.  Min  -4.8621E+00 
Groups 


+1.285E+00 


+5.169E-01 


-1.788E+00 


Figure  D5  Interlaminar  Shear  Stress  (x^j)  in  Segment  of  Warp  Yarn.  Textile  Loaded 
in  Warp  Direction  with  Average  Stress  of  o„°=22.9  GPa 


I 
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Stress  ZZ  (Top) 

Avg.  Max  +8.0458E+00 
Avg.  Min  -i.8510E+00 
Groups 


+6.946E+00 

+5.847E+00 

+4.747E+00 


+3.647E+00 


+2,548E+00 


+1.448E+00 


+3.483E-01 


-7.513E-01 


Figure  D6  Vertical  Stress  (033)  in  Segment  of  Warp  Yarn.  Textile  Loaded  in  Warp 
Direction  with  Average  Stress  of  o,,°=22.9  GPa 


Table  D1  Ratios  of  Maximum  Stress  Components  to  Average  Applied  Warp  Direction 


Component 

Ratio 

_  max/^  0 

011  /011 

2.52 

_  max/^  0 

013  /011 

0.212 

max/  0 

033  /011 

0.350 

REFERENCES 


Pro/MECHANICA,  Release  17.0  (1996),  Parametric  Technology  Corporation,  Waltham 


